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Abstract 


The Finite Element Method (FEM) is extensively used as an engineering analysis 
tool because of its versatility and flexibility. Flowever, the method suffers from 
drawbacks such as discontinuous secondary variables across inter-element boundaries 
and the need for remeshing in large deformation problems. Therefore, researchers in 
recent years have begun to explore the possibility of developing new and innovative 
analysis tools that do not have these drawbacks, and yet have all the advantages of the 
FEM. 

Recent literature shows extensive research work on meshless or element-free 
methods. One such method is the Meshless Local Petrov-Galerkin (MLPG) method. 

This method is based on a local weak form of the governing differential equation and 
allows for a choice of trial and test functions from different spaces. By a judicious choice 
of the test functions, the integrations involved in the weak form can be restricted to 
regular domains. The MLPG method is currently implemented for 2-D potential and 
elasticity problems. 

In this report, the method is further developed for bending of beams - C 1 
problems. A generalized moving least squares (GMLS) interpolation is used to construct 
the trial functions, and spline and power weight functions are used as the test functions. 
The MLPG method for beam problems is applied to problems for which exact solutions 
are available to evaluate its effectiveness. Additionally, a Petrov-Galerkin 
implementation of the method is shown to greatly reduce computational time and effort, 
thus demonstrating that this Petrov-Galerkin approach is preferable over the previously 
developed Galerkin approach. The MLPG method for beam problems yields continuous 
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secondary variables without the need for elaborate post-processing techniques, and the 
accuracy of the method is demonstrated for problems with load discontinuities and 
continuous beam problems. 


This report describes the work that was performed in partial 
satisfaction of the requirements met by Dawn R. Phillips for the 
degree of Master of Science from the George Washington 
University Joint Institute for Advancement of Flight Sciences. 
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Chapter 1: Introduction 


Aerospace structures are very complex in construction. Structural elements used 
are usually built up from doubly-curved shells and stiffeners made up of metallic, 
composite, or sandwich materials. Further, aerospace structures are expected to be 
durable and damage tolerant, are required to have minimum weight, and are expected to 
provide superior performance. These structures are also expected to be in service over a 
wide range of operating conditions and in extreme environments. Satisfying these 
requirements while maintaining cost effectiveness is a complicated but possible task. 

The only efficient way to obtain such a system is through very accurate and high fidelity 
analyses and validation of the resultant design configurations through innovative test 
techniques. 

1.1 Motivation 

The Finite Element Method, because of its versatility and flexibility, is 
extensively used as an engineering analysis tool in civil, automotive, marine, off-shore, 
and aerospace industries. Flowever, the FEM suffers from drawbacks such as 
discontinuous secondary variables (such as stresses) across inter-element boundaries and 
the need for remeshing in large deformation problems. As stresses are discontinuous 
across inter-element boundaries, post-processing techniques are required to achieve 
smooth stress distributions. Four commonly used smoothing techniques are (Cook et al., 
2002) the element smoothing technique, the nodal averaging method, the global 
averaging method, and patch recovery. These methods involve post-processing the FE 
output to obtain smooth secondary variables. 
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The second disadvantage of the FEM is in geometric or material nonlinear 
analysis. In nonlinear analysis, severe mesh distortions can occur. These mesh 
distortions lead to poorly shaped or ill-shaped elements. These ill-shaped elements 
perform poorly and hence remeshing of the deformed analysis region is needed. The 
remeshing and the associated interpolation of the current nonlinear solution onto the new 
mesh is a tedious process. Any method that avoids ill-shaped elements, that provides 
smooth secondary variable distributions, and that retains the advantages of the FEM is 
very attractive. Meshless Methods (MM) appear to show promise in these directions. 

For MM to successfully compete with the FEM, the MM need to be applicable to built-up 
structures. Meshless Methods so far have been applied to one- and two- dimensional C° 
problems. Thus the next step is to apply the MM to C 1 problems involving one 
dimension. In this report, one of the MM, the Meshless Local Petrov-Galerkin ( MLPG) 
method is applied to beam problems. 

1.2 Background 

With the goal of eliminating the disadvantages of the FEM, researchers in recent 
years have begun to explore the possibility of developing new and innovative analyses 
tools that do not have the drawbacks, yet retain most of the advantages of the FEM. 

Nay roles et al. (1992) developed the concept of a diffuse approximation of the 
finite element method. They proposed replacing the finite element interpolation function 
with a smooth function that is diffused and using a moving least squares formulation to 
arrive at the interpolation. A moving least squares (MLS ) interpolation uses the local 
weighted least squares function to evaluate the dependent variable at a point in the 
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domain of the problem. Coefficients in this least squares function are found by 
minimizing the sum of the squares of the error between the interpolation and the value of 
the dependent variable at the nodes. In the FEM, Dirac delta functions are used to 
perform this minimization. For the current Diffuse Element Method, continuous 
weighting functions that vanish at a certain distance from the nodes over which they are 
centered are used. Two very important attributes of the Diffuse Element Method are 
noted by Nayroles et al. (1992): 1) a collection of nodes, without a mesh, and boundary 
conditions are all that are needed to develop the system matrices, and 2) accurate 
solutions are obtained for both regular and irregular nodal spacing. Figure 1 .2.1 shows 
the ways in which domains are modeled in the FEM and MM. 



(a) FEM: nodes and elements 



Figure 1.2.1: Modeling in the FEM and MM 

Belytschko et al. (1994) took the ideas of the Diffuse Element Method further and 
developed the Element Free Galerkin (EFG) Method. In developing their equations, they 
made the important observation that the coefficients in the MLS interpolation should not 
be regarded as constants. As a result, when evaluating the derivatives of the shape 
functions obtained from the MLS interpolation, two very important terms neglected by 
Nayroles et al. (1992 ) were included. The accuracy of the EFG method thus showed 


3 


significant improvement over the accuracy of the Diffuse Element Method. Additionally 
in the implementation of the EFG method, Lagrange multipliers were used to enforce 
essential boundary conditions (EBCs), and a “shadow” cell structure was overlaid on the 
domain to integrate the system matrices. The convergence rate of the EFG method 
depends on the choice of weight function in the interpolation, but significantly exceeds 
that of the finite element method. Several observations were also made about the 
background integration mesh of Belytschko et al. (1994). Because the cells are used 
solely for the purpose of carrying out the numerical integrations, they do not need to 
satisfy the requirements of finite elements, and they can be easily refined in a local region 
(unlike in the FEM). 

Mukheijee and Mukheijee (1997) made important contributions in the imposition 
of essential boundary conditions in meshless methods. They recognized that MLS 
interpolants lack the Kronecker delta property of the usual FEM shape functions. As a 
result, imposition of EBCs is not straightforward. Mukheijee and Mukheijee proposed 
that the values of the dependent variable be replaced by fictitious nodal values to 
accurately satisfy the EBCs at boundary nodes. The resulting system of equations are 
solved for these fictitious nodal values, which are in turn used in conjunction with the 
nodal shape functions to arrive at the numerical solution to the problem. 

While the overlaid cell structure does not have requirements as stringent as the 
finite element mesh, the cell structure is still a mesh that is needed for the EFG models. 
Therefore, one of the advantages of the EFGM is lost. Atluri and Zhu (1998 ) developed a 
truly meshless method that does not require the shadow cell structure to perform the 
numerical integrations. They proposed using a Local Weak Form (LWF), in which 
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calculations begin from the weak form in a local sub-domain. Essential boundary 
conditions are imposed by means of a penalty method. The Petrov-Galerkin method is 
used, as opposed to the Galerkin method used by previous researchers, where the trial and 
test functions are taken from the same space. By a suitable choice of the test function, 
the method can be made local. As such, no overlaying cell structure is required to 
perform the numerical integrations. 

1.3 Objective 

In this report, the MLPG method is first applied to C° problems to understand 
various features of this method. The method is further developed for 1-D C 1 problems 
involving Euler-Bemoulli beams. A Petrov-Galerkin formulation for the beam problems 
is presented. The formulation is applied to several beam problems for which exact 
solutions are available to evaluate its effectiveness. Various features of the method are 
studied and the performance of the method to ranges of important parameters are 
discussed. 

1.4 Scope 

The C 1 problems presented in this report are Euler-Bemoulli beams. Thus, the 
MLPG method is developed using the Euler-Bemoulli beam conventions. These 
conventions are stated as follows: 1) Euler-Bemoulli beams undergo small deformations, 
2) plane sections normal to the neutral axis before deformation remain planar and normal 
to the neutral axis after deformation, and 3) deflection is a function of the axial 
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coordinate alone. A more detailed explanation of Euler-Bemoulli beam theory is 
presented in Chapter 3 . 

1.5 Overview 

In the chapters that follow, the phrase “machine accuracy” appears several times. 
“Machine accuracy” means that the absolute value of the difference between the exact 
solution and the numerical solution is of the order of 10" 14 , using double precision 
arithmetic. 

In Chapter 2, the MLPG method for C° problems, the problems that are described 
by a second order ordinary differential equation, are considered. In C° problems, the 
dependent variables are continuous, but their derivatives may not be continuous. A local 
weak form of the governing differential equation is developed. Approximations to the 
solution known as trial functions are formed using the moving least squares interpolation. 
The Petrov-Galerkin formulation for these C° problems is presented. A system of 
algebraic equations is derived by using the MLS interpolation and the Petrov-Galerkin 
test functions in the local weak form. Numerical examples, including patch test 
problems, mixed boundary value problems, and a typical heat transfer problem are 
worked to evaluate the effectiveness of the method. 

In Chapter 3, the MLPG method for C 1 problems, specifically for Euler-Bemoulli 
beams, is presented. These problems are described by fourth order ordinary differential 
equations. In C 1 problems, the dependent variables and their first derivatives are 
continuous, but higher order derivatives may not be continuous. A local weak form 
(LWF) of the governing differential equation is developed. The moving least squares 
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interpolation scheme is generalized to include derivatives of the dependent variables, and 
is used to construct the trial functions. Test functions are chosen from a different space 
than the trial functions, making the method a Petrov-Galerkin method. The trial and test 
functions are then used in the LWF to derive a system of algebraic equations. 

Numerical examples of beam problems are presented in Chapter 4. A local 
coordinate approach is developed, problem parameters are established, and patch tests are 
performed. Several mixed boundary problems are considered, and the continuity 
requirements for the Petrov-Galerkin test functions are established. Finally, a continuous 
beam problem is studied. 

In Chapter 5, conclusions drawn from the report are presented and summarized. 
Several suggestions for future work are also made. 
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Chapter 2: MLPG for C° Problems 

A Meshless Local Petrov-Galerkin (MLPG) method has been developed for C° 
problems. The method was applied to potential problems by Atluri and Zhu (1998 ) and 
to axisymmetric problems by Raju and Chen (2001). Before C 1 problems can be 
discussed, C° problems must be understood. This chapter presents a description of the 
method applied to C° one-dimensional (l-D) problems. 

First, a local weak form is developed from the classical weighted-residual form of 
the governing differential equation. A moving least squares interpolation is used to 
construct the approximations to the solution known as trial functions. Test functions are 
chosen from a different space than the trial functions, making the method a Petrov- 
Galerkin method. Essential boundary conditions are enforced by a penalty method 
similar to the penalty method employed by the FEM. A system of algebraic equations is 
derived by substituting the trial and test functions into the local weak form. The method 
is evaluated by applying it to several patch test and mixed boundary value problems. 
Finally, a typical example of a heat transfer problem is analyzed using the MLPG 
method. 


2.1 Weak Form for l-D C Problems 


Consider a l-D C problem (Reddy, 1993) governed by 


-J- f Kx) ~ I + c(x)K = fix) 

clx \ clx , 


( 2 . 1 . 1 ) 


in domain Q (0 < x < /) with boundary T, where b(x) and c(x) are problem parameters 
that may be functions of the coordinate x, and /(x) is some “loading,” which may also be 
a function of x. The essential and natural boundary conditions are of the form 
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( 2 . 1 . 2 ) 


u = u on r„ , q = q on T q 

where 


(2.1.3) 

dx 

and T u and T q denote the boundary regions where the primary variable, u, and the 
secondary variable, q, are prescribed, respectively. In 1-D problems, these boundary 
regions are the points x=0 and x=l. The variables u and q represent different physical 
quantities depending on the type of problem considered. For example, in the problem of 
axial deformation of a bar, the primary variable u is longitudinal displacement, b=EA 
where E is the modulus of elasticity and A is the cross-sectional area, /is the applied body 
force on the surface of the bar ( such as friction, self-weight, etc.), and b (du/dx ) , the 
secondary variable, is the axial force. For a heat transfer problem, u is temperature, b is 
the thermal conductivity, /is heat generation, and b (du/dx) is the heat flux (Reddy, 


1993). 


To obtain an approximate solution to Eq. (2.1.1 ), a weighted residual technique is 
employed. As an approximate solution for u is sought, there is an error; that error 
(residual) is 


R = — 


d ( , du\ 


+ cu- f . 


dx v dx ) 

Control of the errors is affected by multiplying the residual by a weight function v(x), 
integrating over the whole domain, and setting the integral to zero: 


(2.1.4) 


f 

d | 

' , dll') ,, 

1 ” 

dx l 

b — +CU-J 

K dx ) 


dx . 


(2.1.5) 
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Equation (2.1.5) represents the classical weighted residual form of the governing 


differential equation. An approximate solution for u is chosen such that each term in the 
approximate solution must be twice differentiable and satisfy all the boundary conditions 
(Eq. 2.1 .2). While these requirements are easy to satisfy in 1 -D problems, for higher 
dimensions, they are difficult to satisfy. Therefore, a formulation that accepts weaker 
requirements on u is sought. The weak form of the weighted residual equation is set up 
by transferring the differentiation from the primary variable u to the weight function v. 
This is achieved by integrating by parts in 1-D and by application of the divergence 
theorem in 2-D and 3-D. Integrating Eq. (2.1.5 ) by parts yields 


0= [b——dx+ [ cvudx- f fvdx- 
J dx dx J J 


Q 


Q 


Q 


, du 
vb — 
dx 


( 2 . 1 . 6 ) 


Integration by parts produces a boundary term [vb- (du/dx )] r . The prescription of the 
secondary variable b-(du/dx) on f is the natural boundary condition (NBC) and is now 
part of the weak form. The requirements on the approximate solution have thus been 
weakened, i.e., u must be differentiable once and must satisfy only the essential boundary 
conditions as the NBCs are included in the weak form. In Eq. (2.1.6), called the weak 
form of the governing differential equation, the chosen approximating functions for u and 
v are called the trial and test functions, respectively. ( The secondaiy variables are 
identified as the coefficients of the weight functions and their derivatives in the boundary 
expressions of the weak form (Reddy, 1993, p. 31).) This weak form is the starting point 
of the Finite Element Method (FEM). 
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In the FEM, u is chosen as a piecewise linear function as shown in Figure 2.1.1. 



Figure 2.1.1: Finite element trial (shape) functions at node j 

The trial functions for element e are chosen as: 

u (e> = N\Ujl\ + N 2 u ( j e) (2.1 .7) 

N el 

where N\ and N 2 are shape functions of the e th element and u = where N e i are the 

e=l 

number of elements in the model. The test function v is chosen as the variation of u: 

v (e) = Su (e) = N\Suj\ + N 2 Su ( j ] . (2.1 .8) 

This choice of v te> as Su {e) makes the FEM a Galerkin method. These choices yield 
several advantages to the FEM: ( 1 ) because the trial functions are piecewise linear, the 
FEM has a local character, and thus the stiffness matrix is banded, (2) the choice v=dn 
yields a symmetric stiffness matrix, and (3) the stiffness matrix becomes positive definite 
after the imposition of boundary conditions because the first integrand in Eq. (2.1.6) 
represents an “energy” quantity. 

The secondary variables are usually the quantities sought in an analysis. For the 
C° problems considered here, the secondary variable is 

q =b— . (2.1.9) 

dx 
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The secondary variable q t for the trial function Uj (see Figure 2.1 .1) is the slope at node /. 
The slopes at node j for elements e and e+1 are obviously unequal. In general, all the 
secondary variables in the FEM are discontinuous across element boundaries because of 
the piecewise nature of the approximation for the shape functions. Post processing 
techniques are required to achieve smooth distributions for the secondary variables. This 
is considered one of the disadvantages of the FEM. 

To overcome the discontinuity problem of the FEM, a diffused element 
formulation was proposed by Nayroles et al. (1992). Later utilizing these concepts, 
Belytschko et al. (1994 ) developed the Element-Free Galerkin method. In these methods 
no elements are present, and trial functions u are formed by passing a smooth function 
through fictitious nodal values (discussed in section 2.2 ). These trial functions are 
written as in the EFG methods as (Mukheijee and Mukhetjee, 1997) 
n 

u(x) = ^ vtj<pj (x) , (2.1.1 Oa) 

7=1 

where n is the number of nodes in the domain of definition of the trial function, U: are 

fictitious nodal values of displacement, and <f>j (x) are shape functions. As the trial 
functions are smooth, the secondary variables are continuous at every point in the domain 
of the trial functions. Using the Galerkin methodology, the test functions are chosen as 
the variation of u, v=Su. and are written in the same manner as the trial functions as 

v(x) = ju\ u) zl U) (x), (2.1.10b) 

where u) u ' are arbitrary constants for displacement, and /^ in are components of the test 
functions. The details of the development of the trial and test functions are discussed in 
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section 2.2. The trial function for node j and test function for node i in the EFG method 
for a 1-D problem are shown in Figure 2. 1 .2. The domain of integration for the i-j term 
in the weak form (Eq. 2.1.6) is the intersection of the trial and test functions and is shown 
by the shaded region, Q. d , in Figure 2.1.2. 


Test function Trial function 



This domain can be large, and its shape may be difficult to determine in 2-D and 3-D 
problems. Because a well-defined shape is desirable for the purpose of integration, a 
background mesh (also called a shadow mesh) - usually rectangular meshes in 2-D 
( Belytschko et al. , 1 994) - is required. As a result, while the formations of the trial and 
test functions do not require elements, the use of a background mesh to perform 
integrations negates the advantage of the EFG method and thus the EFG method is not a 
truly meshless method. 

To develop a truly meshless method, Atluri and Zhu (1998) suggested the choice 
of the test function from a different space, and hence, 

v*Su, (2.1.11) 

and, for example, a weight function whose nonzero values define a well-defined shape 
can be used. Common shapes in 2-D include circles, ellipses, and rectangles. A common 
test function v,- for node i in 1-D (in comparison with a trial function for node j) is 
presented in Figure 2.1.3. These test functions can be chosen to vanish at a certain 
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controllable distance, R 0 , from node i. This localized property of the test functions gives 
the method its local character. 



Figure 2.1.3: Comparison of the domains of the trial and test functions 


Additionally, because the test functions have well-defined shapes and zero value outside 
the local sub-domain £2 S , the integrations can be restricted to £2,. determined from the 
extent of the test functions (see Figure 2.1.3). This choice thus eliminates the need for a 
shadow mesh. The freedom to choose the test function from a different space than the 
trial function makes this a Petrov-Galerkin method. The proposed method is thus called a 
Meshless Local Petrov-Galerkin (MLPG) method (Atluri and Zhu, 1998). 

The weak form is therefore written for the local sub-domain £2 V as 

0= |" b— — dx + \cvudx- \fvdx- vb— 

J dx dx i J |_ dx 

£ 2 , Ll s Ll s 
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As mentioned previously, the test function, v, can be chosen to vanish on r s (see Figure 
2.1 .3). The term [vq]r s is therefore evaluated as zero, and Eq. (2.1.1 6) is reduced to the 
local weak form (LWF) for the MLPG method: 


15 



( 2 . 1 . 17 ) 


f . dn dv , 
b dx + 

icvudx- 

j dx dx 

J J 


£ 2 ,, 


£ 2 , 


£ 2 , 


K] r -\y<i\ 


sq 


where V su represents r v fl T u and F V(/ represents r v H F t/ . The weak form of Eq. (2.1.17) 

is local because the integrations are performed over the local sub-domain Q. s . If the trial 
and test functions of Eq. (2.1 .17) are chosen from the same space via a Galerkin method, 
evaluation of the terms of Eq. (2. 1 . 1 7) yields symmetric stiffness matrices. Thus the 
weak form could be called a local symmetric weak form. (This is the case in the study of 
beam problems by Atluri et al. (1999). ) In this report, a Petrov-Galerkin method is used. 
The resulting stiffness matrices are not symmetric, and thus the term “symmetric” is 
omitted from “local symmetric weak form”. Substitution of the trial and test functions 
into Eq. (2.1.17 ) yields a system of equations of the form 

K ( “ del a +K (Mry) a -f‘“ del =o (2.1.18) 

where the superscript “bdry” denotes boundary, and u are the fictitious nodal values iij 

(see Eq. 2.1.10a). The formation of the system of equations is presented in detail in 
section 2.3. 

Consider now the last two terms of the LWF, 

[v,/] r and \yq ] r . (2.1.19) 

1 su 1 sq 

These terms need to be evaluated at the boundary points. Tire details of these evaluations 
are explained with the aid of a 1-D domain modeled with 17 equally spaced nodes as 
shown in Figure 2.1.4. The nodal spacing in this model is Ax' = //1 6. The primary 
variable, a, is assumed to be prescribed at node 1 , and the secondary variable, q, at node 
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17. In Figure 2.1.4, the test functions are shown at various nodes in the model. These 


test functions are assumed to have an (R 0 / 1) of 2 Ax. 



Figure 2.1.4: Test functions at various nodes in a 17-node model 
Consider the term [vq]r SII - This term must be evaluated for every node in the model 

whose Q s intersects IT. In the model of Figure 2.1.4, there are three such nodes, nodes 1 , 
2, and 3. The key to the contribution of each of nodes 1 , 2, and 3 to the term [vq]r su lies 

in the values of V\, \>i, and v -3 at node 1, where* = 0. First consider node 3: 

V 3 = Oat node 1 , ( 2 . 1 . 20 ) 

and therefore 

Mr,„ = 0. (2-1.21) 

Now consider node 1 : 

Vi = latnodel, (2.1.22) 

and, utilizing Eqs. (2.1.9 and 2.1.10a), 


hdr =Mr = 


, du 
b — 


= b 

d</\ 

d<p 2 

i 

dx 

r 

A X1J 


_ dx 

dx 

dx J 


-lx=0 


U\ 

u 2 


(2.1.23) 
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Finally, consider node 2: 

0 < v 2 < 1 at node 1 , (2.1 .24) 

and therefore 




d<t>2 

dx 


d(t> n 

dx 


x=0 


u 2 


v 2l„0 


(2.1.25) 


Note that the terms h[(d<f>\!dx) (dfo/dx) ... (d<pjdx)]r sii in Eqs. (2.1.23 and 2.1.25) are 

evaluated at node 1 and contribute to the K (bdry) of Eq. (2.1.18) (see Eq. 2.3.1 lb). The 
contribution of node 2 to the term [vq ] r , and ultimately to K (bdry) , is of extreme 

importance and cannot be neglected. 

Now consider the term [va ] r . This term contributes to the f (bdry) of Eq. (2.1.18) 

1 sq 


and must be evaluated for every node in the model whose Q s intersects V q . For a node 
whose v = 0 at node 1 7, 

M r =0. (2.1.26) 

1 sq 

For a node whose v * 0 at node 1 7, [vq ] r is not evaluated as zero unless the 

1 sq 

prescribed secondary variable is zero. The contribution of such nodes to the term 
[vq ] r , and ultimately to the f bdry) , is of extreme importance and cannot be neglected. 

1 sq 

A proper understanding of how the terms of Eq. (2.1.19) are calculated provides 
users of the MLPG method with considerable freedom in choices of nodal spacing and 
sizes of test functions. For the case presented in Figure 2. 1 .4 of a model with equally 
spaced nodes, a choice of a smaller (R 0 / 1) for nodes 2 and N - 1 (for example, here 
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(R 0 / 1) = Ax for nodes 2 and 16 ) ensures that [vq] r = 0 and thus may be preferable. 

However, note that nodes need not be equally spaced. Likewise, the size of Q s for each w 
need not be uniform. When this is the case, a simple choice of a smaller (R 0 / 1) for nodes 
2 and N- 1 may not ensure that all the terms of Eq. (2.1.19) are identically zero for 
additional nodes near the boundaries. In other words, users of the MLPG algorithm 
cannot assume that a simple reassignment of (R 0 /l) will account for the terms of Eq. 
(2.1.19) as in the example above. In order to exploit the full usefulness of the method, 
the terms of Eq. (2.1.19) must be evaluated. 

2.2 Moving Least Squares Interpolation 

Several interpolation schemes are available for constructing trial functions at 
randomly located nodes. The Moving Least Squares (MLS) approximation is one such 
scheme that boasts high accuracy and ease of extension to multi-dimensional problems 
(Nayroles et al., 1993, Belytschko et al, 1994, Atluri and Zhu, 1998, Raju and Chen, 
2001). 

An MLS interpolation is a scheme that passes a smooth function through an 
assumed set of fictitious nodal values. The interpolation is performed such that the least 
squares error between the smooth function and the nodal values is a minimum ( see Figure 

2.2.1) . The MLS interpolations are used to form the trial functions, u, in the current 
implementation of the MLPG method. The trial functions are assumed to be smooth and 
are nonzero over a controllable distance R, from node /'. This distance Rj is usually 
chosen to extend over a much larger extent than the FE shape functions (see Figure 

2.2.2) . The extent of the trial functions can be denoted by Q/j. An MLS approximation 
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can then be made for u\ the value of u in domain Q*. The value for u h is zero outside of 
the domain L>;,. 


U 



Figure 2.2.1: Moving least squares (MLS) interpolation 



The MLS approximation for u in the global domain Q may therefore be written as the 
MLS approximation for u in 12/, as 

m(x) = u h (x) =p T (x)a(x) (2.2.1) 

where 

P T (x) = [pi(x), p 2 (±), p m (x)] (2.2.2) 
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is a complete monomial m th order basis function, and 

a(x) = [aj(x), a 2 (x), •••, a m ( X )] T (2.2.3) 

is a vector of undetermined coefficients. Because the coefficients a(x) may be functions 
of the spatial coordinates, 

x = [x, y, z] T , (2.2.4) 

the values of a(x) can vary with the position of x, y, and z in Q. The global MLS 
approximation is therefore constructed by superposing local MLS approximations in a 
local neighborhood, x of x, where x = x - x . . The local MLS interpolation is then 
written as 

u(x) = Uj(x) = p T (x)a(x) . (2.2.5) 

where p(x) is the basis function, and a( x) and w* (x) are the vector of undetermined 
coefficients and the value of u h (x) in the local neighborhood x , respectively. Examples 
of basis functions for 1-D problems include 

T 

p (x) = [l, x\, linear, m = 2 and (2.2.6a) 


P T (x)= 1, 


X, 


2 

X , 


quadratic, m = 3 . 


(2.2.6b) 


For 2-D problems, basis functions are obtained from Pascal’s triangle (Cook et al., 2002, 
Zienkiewicz and Taylor, 1989) as 


p T (x) = [l, 
p T (x)=k 


X, 

X, 


y]. 


y. 


linear, m = 3 and 

2 2 ~ 
x , xy, y , 


quadratic,/// 


6 . 


(2.2.7a) 

(2.2.7b) 


For 3-D problems, basis functions are obtained from Pascal’s tetrahedron as 
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p'(x) = [], 


X, 


y, z] , linear, m = 4 and 


(2.2.8a) 


T 2 2 2 

p 00= 1 , x, y, z, x , y , z , xy, yz, zx , 


(2.2.8b) 


quadratic,?// = 10 . 

The values of the coefficients a(x) in Eq. (2.2.5) are found by minimizing a 
weighted discrete L 2 error norm defined as (Nayroles et al., 1 992) 


n 

J(x) = ^ Aj (x ) |p ' T (x j )a (Q - u 


7=1 

= [Pa(x)-u ] 1 • ^[Pa(x)-u] 


(2.2.9) 


where A,(x) are weight functions that vanish at a certain distance from x ; -, and n is the 
number of nodes that fall within the local neighborhood x of x,- where Aj ( x ) > 0 . Also 
in Eq. (2.2.9), P is an ( n,m ) matrix, and X is a diagonal (n,n) matrix defined as 


[P] = 


T 

P Ol) 


P 1 (x 2 ) 


T 

P (x») 


( 2 . 2 . 10 ) 


x= 


zij(x) 


Az (x) 


A„(x) 


and 


( 2 . 2 . 11 ) 


U = \u\, U 2 , UnY ■ (2.2.12) 

Note that the values Uj in Eqs. (2.2.9) and (2.2.12) are fictitious nodal values and, in 

general, are not equal to the nodal values of the trial function u h (x) in Eq. (2.2.1) (See 
Figure 2.2.1). 
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Equation (2.2.9) can be written as 


JCO = 


T„T „T 

a P -u 


X ■ [Pa - u] 


= a T P I XPa-2a T P T Xu + u T Xu 


The error norm L 2 is minimized using 
dJ(x) 


daj 


■ = 0 , j = 1,2 . 


Equation (2.2.14) can be rewritten as 
3 J(x) 


da 


= 0 , 


or, 

= 2P T XPa -2P 1 Xu = 0. 
3a 

This leads to 

[A] {a} = [B] {«} 

(m,m)(m,l) (m,ri)(n, 1) 

where 


[A] = P T X P = [B] P =Vi / (T)p(x / )p ( X/ ) 

(m,m) (m,n)(n,n)(n,m) ( m>n )(n,m) “ 

and 

[B] = P 1 X =Ui(x)p(xi), X 2 (x)p(x 2 ), ..., T„(T)p( X// )] . 

(m,n) (m,n)(n,n) 

Solving for {a} in Eq. (2.2.16), 

{a} =[Ar 1 [B] {u} . 

(tn, 1) ( m,m ) (m,n) (w,l) 


(2.2.13) 


(2.2.14) 


(2.2.15a) 


(2.2.15b) 


(2.2.16) 


(2.2.17) 


(2.2.18) 


(2.2.19) 
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Substituting Eq. (2.2.19) into the approximation Eq. (2.2.1 ) yields 


«V) = P T (X>[a] 1 [B] {u} . 

(1 ,m) (m,tn)(m,n)(n, 1) 

The MLS trial functions can then be written as 

n 

u 1 (x) = <I> * (x) • u = ^ <f j(x)u j 
7=1 

where 


( 2 . 2 . 20 ) 


( 2 . 2 . 21 ) 


m 

< I >T (x) = p r ( x)[a] _1 [b] or <t>jW=Yp g W A 'll 

L Jot 


In this report, x = x as 1-D problems are considered. The <pj (x) are called the shape 
functions of the MLS approximation. Also note that <pj(x) = 0 when Aj(x) = 0 (See Eqs. 

2.2. 1 7 and 2.2.1 8). Several weight functions, Aj, were used to construct the trial 
functions, Uj. These weight functions are power weight functions, 


Aj(x) = 



if 0<d;<Rj 
lf d J >R j’ 


(2.2.23) 


where dj is the Euclidean distance between x and Xj denoted by dj = ||x - x 7 -||, and a= 1,2, 
3, and 4, a 3-term spline. 



if 0 <dj<Rj 
if d j > Rj, 


(2.2.24) 


and a 4-term spline. 
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Af(x) = 


1-6 

0 


f d ^ 2 


\ R j J 


+ 8 


< d ^ 

v«jj 


-3 


■Jn 


if 0 < d j < Rj 


if dj>Rj. 


(2.2.25) 


where Rj is a user-defined parameter that controls the extents of the trial functions ( see 
Figure 2.1.3) and is termed the “support of the node j .” (In two dimensions, the “supports 
of the nodal points” are usually chosen as circles of radius Rj.) 

Consider the /V-node model presented in Figure 2.2.3, where N = 9. 



Figure 2.2.3: A 9-node model of a bar 


Figure 2.2.4a presents typical shape functions at nodes.;' =1,3, and 5, evaluated using 
the weight function of Eq. (2.2.23) with a= 4, and Figure 2.2.4b presents the derivative 
d(pj/dx for y=l,3, and 5. These functions were evaluated with a quadratic basis function 
and with (Rj/ 1) chosen as (Rj/ 1) = 0.6. Note that shape functions located equal distances 
on either side of the center nodes of models with uniform nodal spacing are mirror 
images of each other. For example, for the 9-node model presented above, (fh and tfo, fo 
and <pg, etc. are mirror images about the center. 
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2.3 System Equations 

As mentioned previously, the approximations for u are called the trial functions, 
and v are called the test functions. The assumed trial and test functions (Eqs. 2.1.10) are 
substituted into the weak form of Eq. (2. 1 . 1 7), 

0= J b TT dx+ \ cvu dx + a t( M _ U ^ -[v^] r , - \vg ] F ,(2.3.1) 

O s , £2 S , Q s 

to establish the system matrices. The detailed derivation of this system of equations is 
presented below. 

The primary variable, u, is approximated using Eq. (2.2.21): 


l {x) = ^y<pj(x)£j 

7=1 


(2.3.2) 


where are the shape functions, and u , are the fictitious nodal values of u. Substitution 
of Eq. (2.3.2) into Eq. (2.3.1 ) requires the derivative of u h (x). Since u .■ is not dependent 


on x, the derivative is carried out over tpj(x) as 


du y" dipj A 
= > — —U; 

dx dx 

7=1 


(2.3.3) 


The derivative of </>j{x) is obtained as ( Belytschko et. al., 1994) 


tj,x 


B W' + 

8 = 1 


Pg( A 1 b ,x- A lA ,x A ' B )g/J> 


(2.3.4) 


where 



(2.3.5) 
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The steps involved in the evaluation of the derivatives of the shape functions are 
presented in Appendix A. As there are n trial functions used to approximate the primary 
variable, n independent test functions (v,, i = 1,2, ... n) need to be chosen to set up the 
system matrix. Substitution of Eqs. (2.3.3) and (2.3.2) into Eq. (2.3.1) yields 




,u ,■ dx - fvidx 


n, 


7=1 

n 

CA^djitj 
7=1 


fl, 7=1 




(2.3.6) 


-[auvi ] r -MV -MV (/ =1.2, ... n). 

1 su 1 su 1 sq v 


Jr,. 


As discussed in section 2.1, the test functions are chosen as weight functions, 
similar to those presented in Eqs. (2.2.23 - 2.2.25), whose shapes are well-defined. The 
various test functions, Vi, chosen are power functions. 


V/0) = 





i> 


0 


if 0 <dj<R 0 
if dj > R 0 


(2.3.7) 


with dj = |[x — X/|| and / 3 = 1, 2, 3, and 4, a 3-term spline, 


v/(*)H 


1-3 


f d^ 2 


\ R oJ 


+ 2 


A 23 

\ R o J 


if 0 < dj < R 0 


0 


if dj > R 0 , 


and a 4-term spline. 




f dj ^ 

2 

(dA 

3 


(dA 


1-6 


+ 8 



-3 


Vj(x) = < 


k R 0 j 


\ R o J 



V R (l y 


0 


if 0 < dj < R 0 
if dj>R 0 . 


(2.3.8) 


(2.3.9) 
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In Eqs. (2.3.7 - 2.3.9), R 0 is a user-defined parameter that controls the extents of the test 
functions (see Figure 2. 1 .3). A typical plot of the test function of Eq. (2.3.7) with jS= 4 
for node 5 of a 9-node model and ( R 0 /l ) = 2Ax of a bar is shown in Figure 2.3.1. 


1 i 


Zs(?0 0.5 - 


0 



> x 


Figure 2.3.1: Test Function (of Eq. (2.3.7) with J3= 4) at node 5 of a 9-node 
model of a bar 

Substitution of the trial and test functions into Eq. (2.3.6) leads to the resulting system of 
equations 



K(node). +K (b d ry ) ji _f ( no< 1 e ) _ f (bdiy) = Q (2.3.10) 

where the superscript “bdry” denotes boundary, and u are the fictitious nodal values of 
the primary variable u, and 


K (node) 

V 


-f 


, dvj d<t>j 
b- 


Q 


(0 


dx dx 



(2.3.11a) 


and 


„(bdry) 

K ij 



(2.3.11b) 
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and 



S 


and 


(2.3.12a) 


ff hAry) = [auvi ] r (o + [vq ] r (o . (2.3.1 2b) 

su sq 

The “stiffness” matrix K, composed of K (nodc) and K (bdry) , is clearly not symmetric. 
Unsymmetric matrices are not necessarily undesirable. Several numerical methods, for 
example, the boundary element method and the sub-domain collocation method, result in 
unsymmetric matrices. In this meshless formulation, an unsymmetric K is not incorrect 
because, unlike in the FEM, the K matrix in MM is not evaluated from the strain energy 
of the problem, but is obtained by requiring that the weighted residual is zero in an 
integral sense. 

Numerical integration is used to integrate the system of equations as closed-form 
integration of the terms in Eqs. (2.3.1 la and 2.3.12a) is extremely complicated. In the 
Gaussian quadrature integration scheme, an K-point Gaussian will integrate a 2n-\ degree 
polynomial exactly. Equations (2.2.22, 2.3.4, 2.2.17, and 2.2.18) are repeated here for 
convenience: 


<pj(x) = ^pg{ x ) 

g=l 


A *B 


g/ 


(order 2 if quadratic basis is 
used, i.e., if p is quadratic) 


(2.3.13) 


tj,x 


m r 

g = i 


gj+Pgi A l B x -A l A, x A 1 B)^- 


(order 2 ifp is quadratic) 


(2.3.14) 
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(2.3.15) 


n 

[a] = P T IF = [B]P = ^i / (x)p(x / )p T (x / ) 

j = 1 

[b] = P I ) t = U 1 (x)p(x 1 ), ^(x)p(x 2 ), A„(x)iHx n )] (2.3.16) 

The order of Gaussian integration required for acceptable results depends on the basis 
function and weight functions used. The highest order basis function considered is 
quadratic (x 2 ). The highest order weight function available for use as a test function and 
for constructing the trial functions is the weight function of Eq. (2.3.7) with J3= 4, and is 
of the order x 8 . Using this information in Eqs. (2.3.1 la, 2.3.12a, and 2.3.13 -2.3.16), it is 
found that the highest order integrand is of the order x 10 . Therefore, a 6-point or higher 
Gaussian quadrature would successfully integrate the terms of Eqs. (2.3.1 la and 2.3.12a). 
Numerical experimentation showed that an 8-point Gaussian quadrature consistently 
yielded very good results, and is hence used in the numerical implementation of the 
problems presented in section 2.5. 

2.4 Penalty Method for Enforcing Essential Boundary Conditions 

Imposition of essential boundary conditions (EBCs) in the EFG and MLPG 
methods is difficult because the shape functions from the moving least squares 
approximation (discussed in section 2.2) do not have the Kronecker delta property. 
Namely, the Moving Least Squares ( MLS ) shape functions do not pass through the 
fictitious nodal values used to fit them, and unlike in the FEM, 

tj(x k )*8 jk (2.4.1) 

where <p j(Xk) is the shape function for node./ evaluated at nodal point k, and is the 
Kronecker delta. Because the EBCs cannot be directly enforced, a penalty method is 
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employed. In the sections that follow, first, the penalty method in the FEM is explained, 
then, the penalty method used in the MLPG method is presented. 

2.4.1 Penalty Method in the FEM 

In the FEM, a system set of equations is constructed to solve for unknown nodal 
displacements and forces. 

[K]{D}={R} (2.4.2) 

where [K] is the assembled stiffness matrix, {D} is the nodal displacement vector, and 
{R} is the vector of nodal forces. EBCs are input as known displacements, and loading 
and natural boundary conditions (NBCs) are input as known forces. To solve the system 
of equations, the matrices are reordered as 

K uu k un 

_Knu Knn 

where a subscript U denotes values that are unknown, and a subscript N denotes values 
that are known. The resulting equation 

KuuDu + K un Dn=R N (2.4.4) 

can be solved for Du, after which 

K NU D U + K NN D N = R u (2.4.5) 

can be used to evaluate the unknown reactions, Ru. This process of reordering works 
well for small problems and for learning the FEM, but is not used in numerical 
implementation because the process of reordering the matrices requires large amounts of 
memory and run time. A penalty method is therefore employed to solve the system of 
equations. 



(2.4.3) 
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The penalty method in the FEM involves choosing a penalty parameter, or, as a 
very large number (usually 10 20 or 10 30 ). The diagonal stiffness term K u , where i=j 
(corresponding to the known displacement, A), is multiplied by this penalty parameter. 

Similarly, the unknown forces /?, are replaced with aK u bj where ty are the EBCs. This 
inclusion of the EBCs with the force terms rather than with the displacement terms results 
in a system of equations in which the nodal displacements are the quantities sought. 
Consider the / th equation for an M - degree of freedom FE model, 

K n D\ + K i2 D 2 + . . . + K u Di + . . . + K iM D M = R t . (2.4.6) 

This equation can be modified as 

AlA + K i2 D 2 + . . . + aK u Dj + . . . + K iM D M = aK H 5j . (2.4.7) 

The left hand side of Eq. (2.4.7) can be approximated to aK tt Di as this term dominates the 
rest of the terms. Equation (2.4.7 ) can then be written as 

aK ii D l = aKuDj (2.4.8a) 

or 

A =4- (2.4.8b) 

Using this procedure, the prescribed value, 4 > for A is calculated to an accuracy of the 
order (1 la). 

2.4.2 Penalty Method in the MLPG Method 

The penalty method in the MLPG method works in a similar manner to that in the 
FEM. The “assembled” system of equations is 

Ku = f (2.4.9) 

It is desired that 
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(2.4.10) 


a u K u u i = a u K vfii 


or 

ty. u Kjiiii — (x u KjjUj = 0 

a u Ku{ui -U,)=0 ( 2 . 4 . 11 ) 

~ ) = 0 • 

As in the weighted residual sense, because a u (iij -u, ) is not equal to zero, the term is 


multiplied by a weight function v(x) (as in section 2.1) and integrated over the boundary: 



(2.4.12) 


This term for the imposition of the EBCs is included and earned throughout the 
development of the LWF of the governing equation. 

In two- and three-dimensional problems, the boundaries of the domain are 1 -D 
(length) and 2-D (area), respectively, and the integral in Eq. (2.4.12) is evaluated over 
that local boundary segment. In one-dimensional problems, the boundaries are points. 
The integral in Eq. (2.4.12) is evaluated with the dirac delta function as 


ja u (u -u)vS(x = xr m ) dT = [a u (u -t7>] r 


(2.4.13) 


Equation (2.4.13 ) is the form of the penalty method that appears in the development of 
the weak form in section 2.1. Recall the discussion of the terms of Eq. (2.1.19) in section 
2.1. The system of equations is of the form (see Eq. 2.1.18) 


K (nodeA +K (bdry)^_ f (node) _ f (bdry) 


(2.4.14) 
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Also recall the 17-node model of a 1-D domain in Figure 2.1.4, repeated in Figure 2.4.1 


for convenience. 



Figure 2.4.1: Test functions near global boundary 


The primary variable, u, is prescribed at node 1 , where x = 0. Using Eq. (2.1.1 Oa), 

n 

u(x) = ^ <pj ( x)Uj , (2 .4.15) 

7=1 

equation (2.4.13) can be rewritten as 


[a u {u-u)v\ r = a u \ip { fa ... fa\ 

A x Vi 


U\ 

«2 


V Ijc=0 


(2.4.16) 


The term of Eq. (2.4.16) must be evaluated for every node in the model whose Q s 
intersects r„. In the model of Figure 2.4.1, these are nodes 1, 2, and 3. Similarly to the 
terms of Eq. (2.1.19), the key to the contribution of each of nodes 1, 2, and 3 to the term 
of Eq. (2.4.16) lies in the values of V\, v 2 , and v 3 at node 1 . For node 3, v 3 |* = o = 0. For 
node 1, Vi|x = o = 1. For node 2, 0 < V 2 \ x = o < 1. The term of Eq. (2.4.16), evaluated with 
each successive value of v,-| x o for nodes i= 1,2, and 3 contributes to both the K (bdry) and 
the f bdry) of Eq. (2.4.14) (see Eqs. 2.3.11b and 2.3.12b). As previously discussed for the 
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terms of Eq. (2.1 . 1 9), a proper understanding of how the term of Eq. (2.4.1 6) is calculated 
provides users of the MLPG method with considerable freedom in choices of nodal 
spacing and sizes of test functions. 

2.5 Numerical Examples 

In this section, to demonstrate the validity of the MLPG algorithm, the method is 
applied to examples of 1-D C (J problems. The following exact solutions are considered 
for “patch tests”: 

I) u = constant 

II) u — x (2.5.1) 

III) u = x 2 . 

To perform a patch test, each exact solution is prescribed as the essential boundary 
conditions in the problem, and the problem is analyzed with the MLPG algorithm. To 
pass the patch test, the MLPG algorithm must reproduce the exact solution at all interior 
nodes of the model to machine accuracy. In addition to the patch test problems, an 
example problem of heat transfer through a rectangular fin is studied. 

Problem Parameters 

A uniform bar of length / is considered. The bar is modeled using 5, 9, 17, and 33 
equally spaced nodes. The 17-node model is presented in Figure 2.5.1 . 

f* i ►] 

12 9 16 17 

I ^ x 

Figure 2.5.1: A 17-node model of a bar of length / 
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A linear basis function ( 1 , x) should reproduce linear (x 1 ) and lower order 
solutions exactly, and is therefore used for problems I and II of Eq. (2.5.1). Similarly, a 


quadratic basis (1, x, x 2 ) should reproduce quadratic (x 2 ) and lower order solutions 


exactly, and is therefore used for problem III of Eq. (2.5.1). A quadratic basis is also 
used for the heat transfer problem. The weak form (recall Eq. 2.1.17) requires that the 
approximating function, u, be differentiable at least once. The linear basis function is the 
lowest order basis function that meets this requirement, and therefore the lowest order 
basis function that can be used in the MLPG method for C° problems. 

Recall that the governing differential equation is 


d_ 

dx 


^ du' 
b — 
v dx j 


+ cu = f . 


(2.5.2) 


Here, b and c are user-defined constants. The patch tests are performed for various 
chosen values of these constants. 

I. Patch Test - 1: b = 1; c = 0; u = constant = (3\, where jB\ is some arbitrary constant. 
Substitution of these values into Eq. (2.5.2) yields /= 0. EBCs are prescribed at nodes 1 
(x = 0 ) and N (x = P> of an vV-node model as 

”Ui = A 

(2.5.3) 

”U=A- 


This patch test corresponds to an unstressed rigid body displacement (of magnitude (5 \ ) of 
the bar. Values of (R 0 / 1) and (R,-/ 1) were chosen as (R 0 / 1) = 2Ax and (R, / 1) = 1 .0. For 
the 5- and 9-node models, the algorithm calculated the exact solutions for both the 
fictitious nodal values and the interpolated primary and secondary variables. For the 17- 
node model, the algorithm failed to calculate the exact solution for the fictitious nodal 
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values, but the interpolated values were exact. The value of (Rj/ 1) was then reduced to 
(Rj/l) = 8Ax, and with this value and the 17-node model, the algorithm calculated the 
fictitious nodal values exactly. Similar results were obtained for the 3 3 -node model. 
This suggests that the algorithm is capable of reproducing exact interpolated values, but 
exact fictitious nodal values depend on the parameter (Rj/ /). The values of (R 0 / 1) and 
(Rj/ 1) are henceforth chosen as (R 0 / 1) = 2 Ax for all models and (Rj/ 1) = bar length for 
the 5- and 9-node models and (Rj/l) = 8 Ax for the 17- and 33-node models. 

II. Patch Test - II: b = 1; c = 0; u = x/l 

Substitution of b, c, and u into Eq. (2.5.2) yields the loading/ = 0. EBCs are prescribed 
at nodes 1 (x = 0 ) and N (,x = /) of an iV-node model as 


The 5-, 9-, 17, and 33-node models yielded the exact solution with these boundary 
conditions at the nodes and every internal point in the domain, thus passing the patch test. 
The problem can also be worked as the case of a uniform bar with an end load, q (see 
Figure 2.5.2), i.e., with an EBC prescribed at one end and an NBC prescribed at the other 
end. 


/_ 

/ 


~y 


/- 


EA 


A 


Figure 2.5.2: Uniform bar of length / with end load of magnitude*^ 


The prescribed boundary conditions and applied loading are 
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= <7=1 where b = EA 

(2.5.5) 

and 

/ = 0 . 

Once again, the MLPG algorithm with each of the four models reproduced the exact 
solutions for the mixed boundary conditions. 

III. Patch Test - III: b = 0; c = 1 
The exact solution is 



u = 


(2.5.6) 


Substitution of b, c, and Eq. (2.5.6) into Eq. (2.5.2) yields the loading/= (. x/lf . This 
analysis can be performed using three different sets of boundary conditions. 


i) 


ii) 


iii) 


To perform the patch test, EBCs are prescribed at x = 0 and x = I as 


lx=0 


I x=l 


= 0 


= 1 . 


Alternately, mixed boundary conditions are prescribed as 
= 0 


\x=0 


b d U \ 

dx 


= 2bx!l = 0 


X=1 


Thirdly, mixed boundary conditions are prescribed as 
, du I 


dx 


= 0 


I x=l 


x=Q 

= 1 . 


(2.5.7) 


(2.5.8) 


(2.5.9) 
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As expected, the MLPG analysis reproduced the exact solutions for all three cases for all 
nodes of the four models considered. 

Recall the discussions of the boundary terms of Eqs. (2.1.19 and 2.4.16 ). In these 
discussions, it was noted that the size of Q s for each v t need not be uniform and that a 
simple choice of a smaller (R 0 / 1) for nodes 2 and N - 1 may be preferable. For example, 
consider the choice ( R 0 /l ) = 2 Ax for the 1 7-node model of Figure 2.5.1. To account for 
the terms of Eqs. (2.1.19 and 2.4.16), where 0 < V 2 < 1 and 0 < Vi6 < 1, the (R 0 / 1) for 
nodes 2 and 16 is chosen as ( R 0 /I ) = Ax = 0.0625 for a bar of length / = 1 . With this 
choice, the only nodes that contribute to the terms of Eqs. (2.1.19 and 2.4.16) are nodes 1 
and 17. Figure 2.5.3 presents a visualization of the above assignments of (R 0 / /). 




Figure 2.5.3: fl s definitions for various nodes 

The patch tests I, II, and III were performed with these new assignments of (R 0 / 1). As 
expected, the MLPG analysis reproduced the exact solutions to machine accuracy, thus 
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passing the patch tests. These results demonstrate the fact that there is no numerical 
difference between the two choices of ( R 0 /l ), i.e. (R 0 / 1) = 2 Ax uniform for all nodes vs. 
(R 0 / 1) = Ax = 0.0625 for nodes 2 and A-l, as long as the terms of Eqs. (2.1.19 and 
2.4.16) are evaluated correctly. 

In the discussions of the boundary terms of Eqs. (2.1.19 and 2.4.16), it was also 
noted that nodes need not be equally spaced. Consider the 15-node model with unequal 
nodal spacing shown in Figure 2.5.4. 



•• • • • • • 

12 6 12 15 


| > x 

Figure 2.5.4: A 15-node model with unequally spaced nodes 

This model was generated by randomly placing nodes in the region 0 < x < / . The (R 0 / 1) 
for each node was assigned a different value, 

Ax < (R 0 / /) < 2Ax (2.5.10) 

where Ax is the distance between the nodes of the corresponding 1 7-node model with 
equal nodal spacing. For example, for the 17-node model of Figure 2.5.1, Ax = 0.0625. 
The (R 0 / 1) for each node in the model of Figure 2.5.4 was chosen somewhere between 
Ax = 0.0625 and 2Ax = 0. 125. The patch tests I, II, and III were performed with (R 0 / 1) 
for each node assigned as stated above. As expected, the MLPG analysis reproduced the 
exact solutions at all interior nodes in the model and at all interior points in the bar, thus 
passing the patch tests. 
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Example: Heat transfer through rectangular fin 

Consider the rectangular cooling fin shown in Figure 2.5.5. If the variations along 
the _y-direction are negligible, the fin can be modeled as a bar as in Figure 2.5.6, where A 
is the cross-sectional area, P is the perimeter, w is the width, / is the length, and t is the 
thickness. 




x 


Figure 2.5.6: Bar model of rectangular cooling fin 


The governing equation is (Reddy, 1993, pp. 133-134) 

-T^+4(r-r„) = ° 

dx kt 

subjected to boundary conditions 
7'(0) = 7\ Va n 

= 0 


(, A dT 
kA — 
V dx 


X =1 


(2.5.10) 


(2.5.11) 
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where T is temperature, k is thermal conductivity, [3 is the film coefficient is the 

ambient temperature, and Twaii is the temperature of the wall. The equations are recast 
via the non-dimensional quantities 


0 = 


T-T^ 

^Wall ~^oo 



N = 


kt 


(2.5.12) 


as 


^-? + AT0 = O 


di; 

subjected to 

0 ( 0 ) = 1 


(2.5.13) 


r d®' 


yd%j 


= 0 . 


(2.5.14) 


#=1 


The exact solution of the problem is 
cosh N(l - ) 


0 (£) = 0 ( 0 )- 

d® 
dt; 


TT ud® , A r 

H = b— = -bN 


cosh N 1 
sinh N(l - £) 


(2.5.15) 


cosh N l 


In the numerical analysis of the problem, the value of N was chosen as N= 4. The test 
function was chosen as Eq. 2.3.7 with /3= 4. The trial function was constructed from the 
weight function of Eq. 2.2.23 with a= 4 and a quadratic basis function. Tire parameters 
(R 0 / 1) and (R,/l) were chosen as 2 Ax and 8Ax (not exceeding the bar length), 
respectively. The integrations were performed using a 1 0-point Gaussian integration, and 
the penalty parameter was chosen as 10 6 . The bar was analyzed with four models with 5, 
9, 17, and 33 equally spaced nodes. Table 2.5.1 presents the values of the primary and 
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secondary variables obtained with the 5-, 9-, 17-, and 33-node models at stations 
£ = 0, 0.5, and 1.0 along the length of the bar. The values of the exact solutions are also 


included in this table at these stations. All models yielded very good results and the 
accuracy of the solutions improved with model refinement. 


Table 2.5.1: Comparison of the MLPG solution with the exact solution 



Exact solution 
(Eq. 2.5.15) 

MLPG model with: 

5 nodes 

9 nodes 

1 7 nodes 

33 nodes 

m 

<?= o 

1.0 

1.0 

1.0 

1.0 

1.0 

£= 0.5 

0.1378 

0.1360 

0.1377 

0.1377 

0.1378 

<f= 1.0 

0.0366 

0.0420 

0.0369 

0.0360 

0.0366 

de/d £ 

^0 

-3.9973 

-4.1308 

-4.2705 

-4.0322 

-3.9843 

cf = 0.5 

-0.5312 

-0.5502 

-0.5310 

-0.5309 

-0.5305 

<?= 1.0 

0 

0.2737 

0.0468 

-0.0415 

0.0024 


Since the exact solution for this problem is not a simple polynomial, the MLPG method 
did not reproduce the exact solution. Error norms defined as 


and 



I 1 M 

J— ^ (®MLPG ~ ® Exact fg 
\ 2=1 


Ml 


( j M 

II Tf ^ fa MLPG ~ H Exact )g 
1 2=1 


(2.5.16a) 


(2.5.16b) 


were computed at M uniformly spaced points along the bar. These interior points need 
not be coincident with nodes in the model. A value of M= 50 was used. The norms \\e^\ 
and \\e H \\ are presented in Table 2.5.2. As expected, all models yielded accurate solutions 
( within 4%), and the error norms improved with model refinement. 
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Table 2.5.2: Error norm ||g|| for the 5-, 9-, 17-, and 33-node models 


Error norm 

Number of nod 

es in the model 

5 

9 

17 

33 


e& 


0.3127e-2 

0.671 le-3 

0.2195e-3 

0.2154e-4 


e H 


0.3844e-l 

0.6590e-2 

0.5813e-2 

0.3301e-3 


Some post-processing is required to evaluate the secondary variables from the fictitious 
nodal values. To calculate the secondary variables at an interior point, one has a choice 
of two methods. In the first method, the nearest neighboring node to this interior point in 
the domain is evaluated. All the nodes in the domain of influence of this node are 
determined. The nodal shape functions of all these nodes are evaluated at the interior 
point. These shape functions’ values and the fictitious nodal values are then used to find 
the value of the solution u by direct application of Eq. (2.2.21): 

n 

u(x) = ^ Ujtpj(x) . (2.5.17) 

7=1 

Secondary variables may be found in the same direct manner via Eq. (2.3.3): 


du( x ) _ XT' „ dtpj(x) 

dx " J dx 
7=1 


(2.5.18) 


The derivatives of the shape functions are computed at the same time as the shape 
functions themselves, and hence no additional procedures are required. In the second 
method for calculating secondary variables, a shape function is formed over the interior 
point, and all the nodes in the domain of influence of this interior point are determined. 
The fictitious nodal values of these nodes are then used with the value of the shape 
function to find the value of the solution u and the secondary variables via Eqs. (2.5.17 
and 2.5.18). The MLPG and exact secondaiy variable distributions for the 17-node 
model of the heat transfer problem are presented in Figure 2.5.7, and these values agree 
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with the exact solution at all points along the bar. This example demonstrates that one 
can obtain a smooth distribution of the secondary variable. 

k / H 


12 9 16 17 

I >9 



Figure 2.5.7: Comparison of the MLPG and exact secondary variable 
distributions for a 17-node model with uniform nodal spacing 

The same heat transfer problem was then worked using the 15-node model of 
Figure 2.5.4. The ( R 0 / 1) and (Rj / 1) were chosen as in the 1 7-equally spaced nodal 
model. The MLPG and exact secondary variable distributions are presented in Figure 
2.5.8. From this figure, it is seen that the MLPG solution in the region 0<x<//2 is not 
as accurate as the MLPG solution in the region // 2 < x < 1 . This inaccuracy is due to the 
large distance between nodes in the region 0 < x < 1/2 . To improve the accuracy in this 
region, two additional nodes were “sprinkled” into the domain of the problem ( see Figure 
2.5.9). The MLPG solutions before and after model refinement and the exact solution are 
compared in Figure 2.5.9. The inclusion of the two additional nodes significantly 
improves the solution in the region 0 < x < 1 / 2 . 
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Figure 2.5.8: Comparison of the MLPG and exact secondary variable 
distributions for a 15-node model with non-uniform nodal spacing 



Figure 2.5.9: Comparison of the MLPG secondary variable 
distribution before and after model refinement 


2.6 Concluding Remarks 

This chapter presented the MLPG method applied to C° one-dimensional (l-D) 
problems. In the local weak form (LWF) of the governing differential equation, a 
moving least squares (MLS) interpolation was used to form the approximations to the 
solution known as trial functions. Test functions, also needed for the LWF were chosen 
from a different space than the trial functions, making the method a Petrov-Galerkin 
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method. This choice of test functions led to unsymmetric stiffness matrices. The 
essential boundary conditions were enforced by a penalty method, and numerical 
integration was used to evaluate the integrals in the system matrices. The MLPG method 
was applied to and passed several patch test problems. The method was then applied to a 
typical heat transfer problem. Very good results for both the primary and secondary 
variables were obtained. A smooth distribution of the secondary variable was obtained 
without the use of elaborate post processing techniques. 
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Chapter 3: MLPG for C 1 Problems 

In Chapter 2, the MLPG method was studied for the deformation of bars - C° 
problems. In this chapter, the MLPG method is further developed for bending of beams - 
C 1 problems. A local weak form is developed from the classical weighted-residual form 
of the governing differential equation. A generalized moving least squares interpolation 
scheme is used to construct the approximations to the solution known as trial functions. 
Under the Petrov-Galerkin paradigm, the test functions are chosen from a different space 
than the trial functions as combinations of simple weight functions and their derivatives. 
System matrices are derived by substituting the trial and test functions into the local weak 
form. 


3.1 Beam Theory 

The MLPG method for C 1 problems presented in this report was developed using 
the Euler-Bemoulli beam conventions. Consider the beam shown in Figure 3.1.1. Under 
the Euler-Bemoulli bending assumptions, plane sections normal to the neutral axis before 
deformation remain planar and normal to the neutral axis after deformation. The 
deflection w in the z-direction is a function of the x-coordinate alone, i.e., 

w = w(x), u = u (x), and v = 0 . (3.1.1) 


In Figure 3.1.1b, consider AADC in which 



Aw 

Ax 


(3.1.2) 


As Ax -> 0, and for small angles, tan# » 0 gives 



dx 


(3-1.3) 
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where 0 is the slope of the neutral axis. Consider the AABC in Figure 3.1.1c. 
ZB AC = ZCAD = 0 

because normals before deformation remain normal after deformation. 



(a) Beam configuration before and after deformation 




Figure 3.1.1: Euler-Bernoulli beam 


In AABC of Figure 3.1.1c, 


= tan 0 

AB 


or 


(3.1.4) 


(3.1.5 a) 
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u = —z tan 0 ~ — z — — . (3.1.5b) 

t/x 

The strains that correspond to u = -z(dw/dx), v = 0, and w = w(x) can then be evaluated as 


^z = 


dll 

,2 

d w 


dll 

dv 

= 0 

N 

1 

II 

1 ^ 

dx 

Yxy 

dy dx 

()V 



dv 

dw 

■ 0 



II 

h 

dz 

dy 

dw(x) 

0 

y = 

dw 
: h 

du 

dw 

dz 


f zx 

dx 

dz 

dx 


(3.1.6) 


= 0 


Thus all strains except e x are zero. Using the constitutive relationships, the stresses can 
be evaluated. The stress o x corresponding to e x can be evaluated as 


<7 x = Ee x = -Ez d ™ . 

xx 9 

dx 


(3.1.7) 


Now consider the beam segment subjected to a moment in Figure 3.1.2. The moment, M, 
required to return the beam to its undeformed state is 


M 


= ~jV x bdz)z 




Ez 2 b^—ydz 
dx 2 


(3.1.8) 


= E^~ f bz 2 dz 

dx 2 J 
A 


The term J bz 2 dz is the second moment of the area about they-axis and is usually 


termed as the moment of inertia, / 


-yy> 
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(3.1.9) 



z 



Figure 3.1.2: Beam segment subjected to a moment 


3.2 Local Weak Form for Euler-Bernoulli Beam Problems 

The governing equation for an Euler-Bernoulli beam is 
A 

El 4 = / in domain Q ( 0 < x< /) with boundary T (3.2.1) 

dx 

where / is the length and El is the flexural rigidity of the beam, and/is the distributed 
load on the beam. The boundary conditions at x = 0 and x = / can have several 
combinations. The essential boundary conditions (EBCs) are of the form 
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W=w 


( 3 . 2 . 2 ) 


on r w and 


dx 


onTfl 


and the natural boundary conditions (NBCs) are of the form 

V = V on Ty and 

M = M on r M 


(3.2.3) 


where V and M are the shear force and bending moment, respectively, and are related to 
the deflection w as (see Eqs. 3.1.9 and 3.1.10) 


V = -El 


d\v 
dx 3 


and 


M=EI 


d 2 w 

dx 2 


(3.2.4) 


and T w , F s, F v , and F M denote the boundary points where deflection (w), slope ( (9), shear 


( V), and moment (M) are prescribed, respectively. Note that w and V and 9 and M are 
mutually disjoint ( Atluri et al, 1999 and Gu and Liu, 2001 ), i.e., when w = w , the shear 

force Fbecomes the corresponding reaction, and when 0 =9 , the moment M becomes 
the corresponding reaction. 

The weak form of the governing differential equation is obtained in a similar 
manner as for C° problems. The residual error to be minimized is 


R = Eli^-f. (3.2.5) 

dx 

The classical weighted residual form of the governing differential equation for fourth 
order problems is formed by multiplying the residual by a weight function v(x), 
integrating over the whole domain, and setting the integral to zero: 
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v dx . 


(3-2.6) 


0 = 





An approximate solution for w is chosen such that each term in the approximate solution 
must be four times differentiable and satisfy all the boundary conditions (Eqs. 3.2.2 and 
3.2.3). These requirements are difficult to satisfy. Therefore, a formulation that accepts 
weaker requirements on w is sought. The weak form of the weighted residual equation is 
set up by transferring the differentiation from the variable w to the weight function v. 
This is achieved by integrating by parts twice. Integrating by parts once yields 


0 = -EI 


| d ^w dv 

J dx 3 dx 


dx- 



1 

i 

1 / v dx + n x 

EI — -v 

dx 


(3.2.7) 


Q Q 

■7 -1 

where n x [EI(d w/dx )v] p is introduced as a boundary term and n x is the direction 
cosine of the unit outward drawn normal to Q. with respect to the x-axis. The n x thus 
takes values ±1 in 1-D problems. The prescription of the secondary variable 
EI(d 3 w/dx 3 ) on T is a natural boundary condition and is now part of the weak form. 
Integrating by parts a second time to equalize the derivatives of w and v yields 


0 = El P y- v dx- f 
J dx z dx z J 


/ vdx + n x 


El 


d^w 


El 


Jr 


d z w dv 
dx 2 dx 


(3.2.8) 


Jr 


Q Q 

2 2 

where n x [EI(d w/dx )((/v’/T/x)] r is introduced as an additional boundary term. The 

prescription of the secondary variable EI(d 2 v/dx 2 ) on V is also a natural boundary 
condition and is now part of the weak form. The requirements on the approximate 
solution have thus been weakened, i.e., w must now be differentiable twice and must 
satisfy the essential boundary conditions. Additionally, the essential boundary conditions 
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are enforced by a penalty method (Atluri et al., 1 999). As in section 2.4, the penalty 
terms are written as 

a w [(w-w)v] rw (3.2.9a) 

and 


ocq 


f dw ~\dv 


V 


dx 


dx 


(3.2.9b) 


where a w and ae are the penalty parameters to enforce the deflection and slope boundary 
conditions, respectively. Thus, including the penalty terms, Eq. (3.2.8) is written as 


0 = El 




(* d 2 wd 2 v f ^ ] 

I — o -dx- I f vdx + a w [{w-w)v\ Y +a d 

J dx 2 dx 2 J 

a Q. 


dw q ) dv 
dx 


J 


dx 


(3.2.10) 


El 


d 2 w 

dx 2 ’ 


Jr 


El 


d w dv 
dx 2 dx 


In Eq. (3.2.10 ), called the weak form of the governing differential equation, the chosen 
approximations for w are called the trial functions, and v are now called the test 
functions. 

As discussed in Chapter 2, the test functions are chosen independently from the 
trial functions. Test function components chosen in this report for the primary variable w 
in 1-D C 1 problems are the same as those chosen for u in 1-D C° problems. Test function 
components chosen for <9 in 1-D C 1 problems are the first derivatives of the components 
chosen for iv, as 0= dw/dx is also a primary variable (see section 3.3). A typical 
component of the test function v, for node i in 1 -D (in comparison with a trial function 
component (shape function) for node j) is shown in Figure 3.2.1a. As for C° problems, 
these components vanish at a certain controllable distance from node /. 
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5 3 

1^1 —7 ►] 

Domain of the trial function (2i? ; ) 

(a) Deflection, w 


Component of 
test function 



(b) Slope, 6 

Figure 3.2.1: Comparison of the domains of the trial and test functions 


The derivatives of these components also vanish at the same distance from node i (see 
Figure 3.2.1b). This localized property of the test functions preserves the local character 
of the method. The integrations over £1 become integrations over a local sub-domain, O v , 
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and the Q v can be determined from the extent of the test functions (see Figures 3.2.1). 
The weak form is therefore written for the local sub-domain Q v as 



where r w and T s e are the boundaries where w and 0 are prescribed on the local boundary 
f T s fl T w and Y s D Y s ). Note that if the local boundary Q v does not intersect F w or F e 


(i.e. when the Q. s is completely within the interior of Q), the penalty terms are not 
considered for that local boundary. Recalling Eqs. (3.2.4), Eq. (3.2.1 1) is written as 



When the local boundary T s intersects the global boundary F. four boundary conditions 
are possible ( Atluri et ah, 1999): 

r,nr„,, r,nr„. 

(3.2.13) 

r s . nr F , and r, n y m . 

Utilizing these subsets, Eq. (3.2.12) becomes 


57 




As mentioned previously, the test function, v, and its derivatives can be chosen to vanish 
on T s (see Figures 3.2.1). Equation (3.2.14) then is reduced to the local weak form 
(LWF ) for the MLPG method: 



where, as in Eq. (3 .2 . 1 1 ), r. w represents T s f| r u , and Y s e represents F v f| Tg , and 
similarly, IV represents r v fir F and TW represents r v f| T v/ . Now n x is the direction 
cosine of the unit outward drawn normal to Q. s ; n x = 1 if the boundary is on the right side 
of Q s , and n x = -1 if the boundary is on the left side of Q. s . The weak form of Eq. (3 .2. 1 5 ) 
is local because the integrations are performed over the local sub-domain Q s . 

The trial functions are written as 

n 

w(x) = (w / yr ( / M,) (x) + OjWp (*) ] , (3.2.16a) 

7=1 
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and the test functions are written as 


v(x) = Pi W) Zi W) (x) + pj 0) Zi 0 hx) ■ 


(3.2.16b) 


As discussed in Chapter 2, if the trial and test functions of Eqs. (3.2.16) are chosen from 
the same space via a Galerkin method, symmetric stiffness matrices are obtained from 
Eq. (3.2.15). Again, this is the case in the study ofbeam problems by Atluri et al. (1999). 
In this report, a Petrov-Galerkin method is used, and thus the resulting stiffness matrices 
are not symmetric. The details of the development of the trial and test functions are 
presented in sections 3.3 and 3.4. Substitution of the trial and test functions into Eq. 
(3.2.15) yields a system of equations of the form 

K°” de) d + K (bdry, d-f (nodel -f (bdry) =0 (3.2.17) 

where the superscript “bdry” denotes boundary. Note that the locality of the MLPG 
method (as integrations are performed over £2,) makes the stiffness matrices of Eq. 
(3.2.17) banded. This is one of the advantages of the FEM that is retained by the MLPG 
method. The detailed formation of the system of equations of Eq. (3.2.17) is presented in 
section 3.4. 

3.2.1 Boundary Terms in the LWF 

As in Chapter 2, the boundary terms in the weak form need special attention. The 
issues related to these boundary terms are discussed below. 

Consider the boundary terms of the LWF: 


,3 

T7T Cl w 

El — — v 
dx 


El 


d 2 w dv 


dx 


2 dx 


l s8 


(3.2.18a) 
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(3.2.18b) 


, n r 

M— 

sV ’ x 

dx 

w)v] r 

A SW 

, ae 


L sM 


dw q ' dv 
dx J dx 


Jr, 


s(9 


(3.2.18c) 


The terms of Eq. (3.2.18a) resemble the term [v</]r iH of Eq. (2.1.19), and the terms of Eq. 
(3.2.1 8b ) resemble the term \\ q ] r of Eq. (2. 1 . 1 9). Likewise, the terms of Eq. (3 .2. 1 8c) 


resemble the terms of Eq. (2.4.1 6). These terms need to be evaluated at the boundary 
points. The boundary term evaluations are explained with the aid of a typical 1 7-node 
model of a beam as shown in Figure 3.2.2. 



The primary variables, w and 6, are assumed to be prescribed at node 1 , and the 

secondary variables, V and M, at node 1 7 . Recall that w and V and 0 and A/I are 
mutually disjoint, i.e., for example, w and V cannot be prescribed on the same boundary 
point. In Figure 3.2.2, the test functions with an (R 0 /l) of 2Ax located at nodes 1, 2, 3, 

and 9 in the model are shown. Consider the term n x [eI [d '\\ j dx ' )rj r of Eq. (3 .2. 18a). 

This term must be evaluated for every node in the model whose Q. s intersects T w . In the 
model of Figure 3 .2.2, there are three such nodes, nodes 1 , 2, and 3 . The key to the 
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contribution of each of nodes 1,2, and 3 to the term n x El [d ’ vr/ dx' )v r lies in the 


values of vj, v 2 , and v 3 at node 1 , where x = 0 and n x = -1 . First consider node 3: 


v 3 = 0 at node 1 , 
and therefore. 


(3.2.19) 


El 


,3 

a w 


dx 


-v 3 


= 0 . 


Jr„ 


(3.2.20) 


Now consider node 1 : 


V\ = 1 at node 1, 

and, substituting Eq. (3.2.16a) into the term n x \El(d J w/ dx ' )v 


(3.2.21) 


El 


d^w 
dx 


■Vi 


= n , 


J r s , 


El 


,3 

d w 
dx' 


-1EI 


,3 (w) 3 (w) 

d y/\ dy/2 


dx 


dx 


-1 -El 


,3 (9) ,3 (9) 

d y/| d y/ 2 


dx 


dx' 


,3 (w) 
d ¥n 


dx 


,3 (<9) 
^ Vn 

dx 


-h=0 


x=0 


H’j 

w 2 

#2 


(3.2.22) 


Finally, consider node 2: 

0 < v 2 < 1 at nodel, (3.2.23 ) 

and therefore, 
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(3.2.24) 


El 


d 3 w 


dx 



d^Y { w) d ? 'i// ( 2 w) d\j/ { n 


= -l-EI\ 


dx 


dx 


\,A W ) 


dx 


3 


x=0 


W\ 

w 2 

w„ 


n3*3 


x=0 


+ 3'xT 


x=0 


■l-EI 


dV x 0) d^Y^ 


dx 


dx 


3 


dYP 


dx 


3 


Jjc=0 


Ox 

&2 

o„ 


333 


x=0 


+4PAP 


x=0 


Note that the terms 


-\-EI 


3 (w) ,3 (w) 

d y l dYi 


dx 


dx" 


3 (w) 
d Yn 

dx 


-U=0 


and 


-\-EI\ 


3 (0) 
d Y\ 

dx 


3 (0) 
d Y2 

dx^ 


3 (0) 
d Y n 

dx 


Jjc=0 


in Eqs. (3.2.22 and 3.2.24) are evaluated at node 1 and contribute to the K (bdry) of Eq. 
(3.2.17) (see Eq. 3.5.4e). 

The remaining terms of Eqs. (3.2.18) are evaluated in the same manner as the 
terms of Eqs. (2.1.19 and 2.4.16) and using the trial and test functions of Eqs. (3.2.16 ) as 

discussed above. Consider the term n x \El[d 2 w/dx 2 \dv/dx)] r g of Eq. (3.2.18a). This 


term must be evaluated for every node in the model whose Q. s intersects T e . In the model 
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of Figure 3.2.2, these are nodes 1 , 2, and 3. The key to the contribution of each of nodes 
1, 2, and 3 to the term h x \eI { d 2 w/dx 2 )(rfv/c&c)J r lies in the values of (dvddx), (dv 2 /dx), 


and (dvddx), at node 1, where x = 0 and n x = -1 . First consider node 3: 
(dvddx) = 0 at node 1 , 
and therefore, 


EI\ d \\j dx 


(dv 3 /dx) 


= 0 . 


1 sff 


Now consider node 1 : 

Vi = 1 and (dvddx) = 0 at node 1, 
and therefore, 


EI\ d 2 w/dx 2 


{dvi/dx) 


= 0 . 


Jr, 


s6 


Finally, consider node 2: 

(dv 2 /dx) is nonzero at node 1 (in fact, dv 2 /dx < 0 in Figure 3.2.2) 
and, substituting Eqs. (3.2.16) into the term n x [El[d 2 w/dx 2 \dv/dx) 




(3.2.25) 


(3.2.26) 


(3.2.27) 


(3.2.28) 


(3.2.29) 


63 



(3.2.30) 


= -l-EI 


~ d 2 ¥ [ w) dV 2 w) 


dVn W) 


Wo ..(w) d%2 
. dx 


,(d) d%2 


1 El\ d2y/({0) 


d 2 Vn 6) k \ u {w) d X2 

,2 ^2 r. 


Ad) d X\ 


Note that the terms 




,2 (w) ,2 ( w ) ,2 (w) 

dy/\ dy / 2 d y/ n 


,2 {0) ,2 (d) ,2 ( 6 >) 

dy/\ dy /2 d yr n 


in Eq. (3.2.30) are evaluated at node 1 and contribute to the K y of Eq. (3.2.17) (see 
Eq. 3.5.4e). 

Now consider the term n x lV l ’Jr vF of Eq. (3 .2 . 1 8b ). This term must be evaluated 

for every node in the model whose Q. s intersects T v , where x = / and n x = 1 . For a node 
whose v = 0 at node 1 7, 


: [v v Jr vl/ 


( 3 . 2 . 31 ) 
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For a node whose v ^ 0 at node 1 7, n x [v v lr> is not evaluated as zero unless the 
prescribed shear is zero. Substitution of Eq. (3.2.16b) into the term n x IV 1 Jr . yields 


VvL =v f - (w) - (w) 

If vV 


Mi Xi 


x=l 


( 9 ) ( 6 ) 

+ Mi Zi 


x=l 


(3.2.32) 


Similarly, the term n x Midv/dx )]p of Eq. (3.2.18b) must be evaluated for 

every node in the model whose Q. s intersects V M , where x = / and n x = 1 . For a node 
whose (dv/dx) = 0 at node 17, 

n x [M(dv I dx)\- sM =0. (3.2.33) 

For a node whose (dv / dx) * 0 at node 17, n x \M ( dv / dx ) r is not evaluated as zero 

L JA sM 

unless the prescribed moment is zero. Substitution of Eq. (3.2.16b) into the term 
n\M(dv I dx)\- yields 


Midvldx) L = M 

1 sM 


Mi 


(V1-) dxl 


(w) 


dx 


+ Mi 


(8)dZi 




x=l 


dx 


x=l J 


(3.2.34) 


Note that the terms 


Xi 


(w) 


X—l 


Xi 


m 


dX\ 


(w) 


x=l 


dx 


,and 


dXi 


(0) 


dx 


x=l 


x=l 


in Eqs. (3.2.32 and 3.2.34) are evaluated at node 17 and contribute to the f ,bdry) of Eq. 
(3.2.17) (see Eq. 3.5.4g). 

Now consider the penalty term a w [(vi’ - w)v] r of Eq. (3.2.1 8c). This term must 


be evaluated for every node in the model whose Q. s intersects r w . Again, these are nodes 
1 , 2, and 3. The key to the contribution of each of nodes 1 , 2, and 3 to the term 
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a w \{w - w)v] r lies in the values of V\, v 2 , and v 3 at node 1, where x = 0. Substitution of 

Eq. (3.2.16a) into the penalty term yields 
a w [{w-w)v] rw 


= a 


w 


( w ) (w) (w) 

W\ V 2 ■■■ Vn 


M’| 

W 2 


V, 


Jx=0 


d x =o 


w„ 


(3.2.35) 


+ O' 


(0) (d) ( e ) 

¥\ V 2 ••• ¥„ 


Jx=0 


*1 
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V, 


dx=o 


Or, 


■ a w w- v, 


dx=0 ' 


For node 3, v 3 |a: = o = 0. Fornode l,Vi| x = o = 1. The term of Eq. (3.2.35) is evaluated with 
each of these values. Fornode 2, 0 < v 2 |x = o < 1, and substitution of Eq. (3.2.16b) into Eq. 
(3.2.35) yields 
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(3.2.36) 


a w [(w-iv)v] r 

L SW 










This term must be evaluated for every node in the model whose C2 S intersects F g. Again, 
these are nodes 1, 2, and 3. The key to the contribution of each of nodes 1 , 2, and 3 to the 
term a 0 i(dw/ dx)-Q \dv/ dx )j r ^ lies in the values of (dvddx), (dv 2 /dx). and (dvddx), at 
node 1 , where x = 0. Substitution of Eq. (3 .2. 1 6a) into the penalty term yields 
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For nodes 3 and l, [dvt/dx] x = 0 = 0(i = 1,3), and therefore 

ocq {{ K dw/dx)-Q\dv/dx)\ r ^ = 0 . For node 2, (dv 2 /dx) is nonzero, and substitution of Eq. 
(3.2.16b) into Eq. (3.2.37) yields 
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( dw ~ dv 

a e 0 — 

v dx ) dx 


dy/\ w) dy/ ( 2 w) di// ( n w) 


<u y 2 


, ( w ) 
(w) d% 2 


(9) d% 2 


, «9) , (#) 

dy/ x dy/ 2 


, (0) 
dVn 


B l ( A ^ W) 

a ( w)dl 2 
* 2 ^ ~ 


W)dx 2 
M2 «fr 




, (9) 

. . .(g) d X2 

2 dx 


The teims 


d y/\ w * d \f/2 * 


. (9) , (9) , (9) 

d y/\ d y / 2 d y/ n 

dx dx dx 


of Eqs. (3.2.35 - 3.2.38) contribute to the K (b y) of Eq. (3.2.17) (see Eq. 3.5.4e). The 
terms 


(H-) (9) dx) 

/Ci > /Ci * , 

X=l x—l dx 


of Eqs. (3.2.35 - 3.2.38) contribute to the f ( •’ of Eq. (3.2.17) (see Eq. 3.5.4g). 
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As discussed in Chapter 2, a proper understanding of how the terms of Eqs. 
(3.2.18) are calculated provides users of the MLPG method with considerable freedom in 
choices of nodal spacing and sizes of test functions. The value of (R 0 / 1) maybe adjusted 
in certain cases to account for the terms; however, in order to exploit the full usefulness 
of the method, the terms of Eq. (3.2.18) must be evaluated. 

3.3 Generalized Moving Least Squares Interpolation 

Recall from section 2.3 the MLS interpolation scheme for constructing trial 
functions for C° problems. The local MLS interpolation is written as 

’T 

u(x) = «x( x ) = P (x)a(x) where p(x) is the basis function, and a(x ) is the vector of 
undetermined coefficients in the local neighborhood x . The values of a(x) are found by 

minimizing a weighted discrete L 2 error norm. The 1-D shape functions resulting from 

n 

this MLS interpolation scheme are u(x) = £0j(x)Uj . Note that only the interpolation 

7=1 

for the primary variable, displacement, is performed. 

In beam problems, both the deflection w and the slope 6 are the primary variables. 
In the FEM, the Hermite functions are used as interpolation functions for the primary 
variables. See Figure 3.3.1 for a comparison of the FEM Lagrangian and Hennite shape 
functions. The additional information (i.e., the slope) used in the Hennite shape 
functions must also be used in the approximation of the MLPG method. In order to 
accomplish this, a generalized moving least squares (GMLS) approximation is developed. 
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U — (f)\ t/j + (J>2 U 2 


W = 0 l W 1 + <p 2 W 2 + 0 3 H’ 3 + (f> A W A 


W\ H’2 



(b) C 1 problems - Hermitian functions 


Figure 3.3.1: Comparison of FEM shape functions for C° and C 1 Problems 

In this report, the spatial coordinates y and z of x = [x y z] 1 are not present as 
1-D problems are considered, and therefore x =x. The GMLS approximation for w in a 
global domain il may be written as in the MLS procedure as 

w(x) = w 1 ( x ) = p T (x)a(x) . (3.3.1 a) 

Likewise, the local GMLS approximation is written as 

Mx) = wj(x) = p T (x)a(x) (3 3 . lb) 

where x = x-Xj , p(x) is the basis function, and a(x) , the vector of undetermined 

coefficients, is found by minimizing a weighted discrete H h error norm (Nayroles et al., 
1992, Atluri etal. , 1999): 
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H h (a) = ^ ^ Aj (x)[D a p T (xj )a -D 
7=1 


H’,- 


(3.3.2) 


where A/ is a weight function, D" denotes the 0 th derivative, and X indicates the 

\a\<h 

summation of all derivatives up to order h. 

For beam problems, the primary variables are the deflection w and slope 9 = 
(dw/dx), and hence, the weighted discrete H 1 error norm is used: 


H l (a) = V V ^ (x)[D a p T (xj )a - D a wj 
7=1 |ar|<:i 


-2 

7=1 


/y" ! (.v) p 1 (.Vy )a - Wj 


+ A { f ) (x) 


d]i'(Xj) 

dx 


a-9j 


(3.3.3) 


In this report, Atp (x) and X ( f\x) are chosen as the same weight functions, and will 
hereafter be referred to as A Ax ) . In matrix form. Eq. (3.3.3) is 


H [ ( a) = [Pa - w] 1 A[Pa - w]+ [p r a -t] 1 ^[P A a -t] 

fA-ilTr, 


W 


A 0 

0 A 


w 


(3.3.4) 


= [Qa - s] T A[Qa - s] 

where P and P Y are ( n,rn ) matrices and A, is a diagonal (n,n) matrix defined as 


[P] = 


IT T 

P Ol) P ( x 2) ••• P ( x n) 


(3.3.5a) 


[Pj = 


XT X 

PxU'l) P.r( x 2> ••• Px( x n) 


(3.3.5b) 
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M*) 


h(x) 


Xyj (A) 


where 


(3.3.6) 


T 2 m - 1 

p (x) = 1, X, X , ... X 


T c/p (x) [ m-2 

Px (x) = — = 0, 1, 2x, ... (m-l)x 

dx 


with (m- 1) as the order of the 1-D basis function p(x) used in the MLS approximation. 


I 0 
0 k 


are the basis function matrix, the nodal displacement vector, and the weight function 
matrix, respectively. Further manipulation of Eq. (3.3.4) leads to 

H l =[a T Q T -s T jA[Qa-s] 

= a 1 Q t A-s 1 A][Qa -s] 

(3. 

= a 1 Q r AQa -a r Q 1 As -s 1 AQa +s T As 
= a 1 Q r AQa -2a T Q T As + s T As. 


The norm H 1 can be minimized using: 


= 0 , / = 1 , 2 , ... « . 


Equation (3.3.10) can be rewritten as 
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(3.3.11) 


9a T 

9/7 ^ T t ~ 

= 2Q i AQa-2Q 1 As = 0. 

9a T 

Equation (3.3.1 1) leads to 

[A] {a}= [B] {s} 

(m,m)(m,\) ( m,2n)(2n,\ ) 

where 

[A] = Q T A Q 

(m,m) (m,2n)(2n,2n)(2n,m) 

= P T X P + pj X P Y 

(m,n) (n,n) ( n,m ) ( m,n ) (nji) ( n,m ) 

and 


(3.3.12) 


(3.3.13) 


[B] = q’ 1 A = 

( m,2n ) (m,2n) (2n,2n) 

Solving for {a} using Eq. (3.3.12) gives 
{a} = [a] -1 [B] {s} 

(m, 1) (m,m) (m,2ri)(2n,\) 

Substituting into the approximation Eq. (3.3.1a) 
w / 7 (x) = p T (x)[a ] _1 [b] {s} . 

(1 ,m) ( m,m ) (m,2ri) (2n,l ) 


p r X , P Y X 

(m,n)(n,n) ( W , H )(w,w) 


(3.3.14) 


(3.2.15) 


(3.3.16) 


The trial functions used for beam problems are finally written as a linear combination of 
nodal shape functions: 


w(x) = [wjy/f ) (x) + 0jvf\x)j 

7=1 


(3.3.17) 
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where 


V 


r ) M=£vv[ A_lpT *-]® 


g = 1 


(3.3.18) 


jfl 

vf\x) = ^ '[ A “ lp J >-],/■ 


g = 1 


Note that w.- and Q j in Eq. (3.3.17) are fictitious nodal values of deflection and slope, 
respectively. 

As in Chapter 2, three types of weight functions X,(x) are considered for 
constructing the trial functions: power functions, 

i 


l 7 -(x) = 


a 


if 0 <dj<Rj 
if dj>Rj, 


(3.3.19) 


with dj = \\x -Xy||, the Euclidean distance between x and Xj. and a= 1 , 2, 3, and 4, 


/Ux) = 


1-3 

0 


f d> 2 

Jn 


+ 2 


\ T i) 


if 0 < dj < Rj 


if d j > R j. 


(3.3.20) 


a 3-term spline with dj=\\x- Xj ||, and a 4-term spline. 


A/(x) = 


1-6 

0 




\ R J J 


+ 8 


■J.n 


-3 


<d, ^ 


if 0 < dj < Rj 


if d j > Rj 


(3.3.21) 


where Rj is a user-defined parameter that controls the extents of the trial functions. These 
weight functions are chosen to demonstrate the robustness of the MLPG method. 
Consider the 17-node model presented in Figure 3.3.2. 
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5 


9 


16 17 


1 2 

| > x 

Figure 3.3.2: A 17-node model of a beam of length 4/ 



0 x/4 1 1 


(a) Shape functions, y/ 9 



0 x/4/ 1 

(b) Derivatives of the shape functions, dyr 9 / dx 

Figure 3.3.3: Typical shape functions and their derivatives 

Figure 3.3.3a presents y/ ( p and y/f ] at node j = 9, typical shape functions evaluated 
using the weight function of Eq. (3.3.19) with a= 3. Figure 3.3.3b presents the 
{dy/p / dx) and the (dy/'p ' dx) for node / = 9. These functions were evaluated using a 
quartic basis function. The value of (R,/ 1) was chosen as (Rj/ 1) = 3.5. From these plots, 
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it is seen that while 6 = dw/dx, y/ H ' (dy/ [ w) / dx ) . This is because the derivative of 
y} w) involves the inverse of the [A] matrix. One should note that, while 
y/ u>) -£■ ( dy/ (w) /dx) , the basis function used for y} 0) must be the derivative of the basis 
function used for y} w> . For example, if a quadratic basis (1 , x, x 2 ) is used for y)"\ then (0, 
1 , 2x) must be used for yi H ' . These are important characteristics of the MLPG method for 
Euler-Bemoulli beam problems. 

3.4 Test Functions Used 

The MLPG equations are derived using the weak form of the governing equation. 
Recall the weak form from section 3.2: 



Also recall Eq. (3.3.17), in which the trial function w is approximated as 

M.x) = ^(wjvy\x) + 0jysf\x)) (3.4.2) 

7=1 

where Wj and 6j are the fictitious nodal deflections and slopes of the trial function, 

and yi w) and yi 0] are their corresponding shape functions, respectively, given by Eqs. 
(3.3.18). The test function, v, is approximated using 

v(x) = (^) (3.4.3) 
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where //• " ) and //■ 6 are the arbitrary constants for deflections and slopes of the test 
function, and xP( x ) and xf\ x ) are components of the the Petrov-Galerkin test 
functions that are chosen from a different space than y/-' ] (x) and y/f ) (x) . Recall that 
the expressions for the shape functions are written as in Eq. (3.3.1 8) as 

m r , 


V 


(w> (x) = ^pg(xj 




gf 


g= 1 


(3.4.4) 


V 


( 0 ) 


m 


(*>=!>* 

g= 1 


(Xj) 


- l pjl 


gf 


In a Galerkin approximation, the components of the test functions Xi ( x ) an ^ Xi ( x ) 
would take on the same form as the shape functions: 


m 


x { f ){x ) = Yu p g {x j 

g = i 


m 


xf\ x ) = ^p g ( X j 

g = i 


|A -1 P -1 ' X 


|A _1 PJ 


£7 




(3.4.5) 


As mentioned previously, Atluri et al. (1999) used the Galerkin approximation of Eq. 
(3.4.5 ). However, in this work, a Petrov-Galerkin approximation is used and the 

components of the test functions %\ w) (x) and %) 0) (x) are chosen from a different space 
than the shape functions y/J :) (x) and y/ ( p(x) . The Petrov-Galerkin formulation is 

further discussed in section 3.6. The test function components x\ w) (x) are chosen as 
simple weight functions similar to those of Eqs. (3.3.19 - 3.3.21 ) as 
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zrM= 


\-[d}l R l) 


\p 


if 0 < dj < R 0 


if dj > R 0 


(3.4.6) 


with di = |[x -x,-|| and 0=2,3, and 4, 


*i w, OOH 


f ,7 A 2 


1-3 


0 


dj 

\ R oJ 


+ 2 


' d X 


\ R o y 


if 0 <dj<R 0 


if dj > R 0 


with dj = | [x — X;||, and 


( 3 . 4 . 7 ) 




1-6 

0 


'd,' 2 


\ R 0 J 


+ 8 




\ R 0 J 


f 1 \ 


-3 


d 


\ R 0 J 


if 0 < dj < R 0 


if dj>R 0 . 


(3.4.8) 


where R 0 is a user-defined parameter. Plots of the components of the test functions 
x\ w \x ) and Xi &> ( x ) of Eq. (3.4.6) with 0=4 for node 9 of a 17-node model of a beam 
with ( R 0 /l ) = 2A x are shown in Figure 3.4.1a. The corresponding derivatives, 

(dx ! M * / dx) and (dXj d) / dx) are shown in Figure 3.4.1b. Note that for the test 


functions, xf^ = (d%\^ I dx) . 
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Figure 3.4.1: Typical test function components and their derivatives 

3.5 Development of the MLPG Equations 

To evaluate the integrands and the terms involved in the weak form, the 

derivatives of w and v, the trial and test functions, are needed. Since w f , 6 j , jU - W) , and 
ji\ B) are constant values, the derivatives are earned out over y) w) , y) 0> , ^ w \ and ^ H) as 
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and 


dw _ y 1 

//r _ 


dx 


w 

dx 2 


7=1 


J dx J dx 
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n 

=E 


f a.>o 


,3 » 

d w 


dx 


7=1 

n 

■I 

7=1 


w 


<7 XI/, 

J + 0 , 


,2 (0)^ 

d y/ j 


J dx 2 


J dx 2 


w 


(P Xj/^ , 


,3 (0) 
d y/ 


J dx 


] dx } 


, , (w) J (0) 

= (w) dx t (0) dXi 
dx * dx ' dx 

,2 ,2 ( w ) ,2 ( 0 ) 

d v _ j w ) d Xi , „«?) d Xi 
2 Ai 2 2 

dx rfx <7x 

Appendix A presents explicit expressions for all the derivatives of y / w) and 
Substitution of Eqs. (3.3.17, 3.4.3, and 3.5.1) into the weak form leads to 


(3.5.1a) 


(3.5.1b) 
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(3.5.2) 
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Requiring that Eq. (3.5.2) be valid for arbitrary values of // ; (in and /// <9) leads to the 


MLPG equations as 

K (node)^ +K (bdry)^ _ f (node) _ f (bdry) _ 
where “bdry” denotes boundary and 


(3.5.3) 
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(3.5.4a) 


d = {w, 9 X , w 2 , 9 , ...} 

are the fictitious nodal values of deflections w and slopes 9, and 


with 
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(node) 
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(bdry) 


(node) 
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kf de) =£/ 
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d 2 x m 


dx 
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-dx 
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dx 


dx 
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dx 
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dx 


dx 


(3.5.4b) 

(3.5.4c) 

(3.5.4d) 
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(3.5.4e) 


(3.5.4f) 


(3.5.4g) 
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where i,j =1,2 , n. 


After these equations are solved for the fictitious nodal values Wj and 9 j , the 

interpolated primary and secondary variables may be computed. The deflection, w, at 
any point in the beam is calculated from Eq. (3.3.17), 

n , 

w(x) = y ] [wjYj (x) + OjY] (x ) j , (3.5.5) 

7=1 

and post processing is accomplished by either of the methods discussed in section 2.5. 
The slope 6 , moment M, and shear V can just as easily be calculated from Eqs. (3.5.1a). 

3.6 The Petrov-Galerkin Formulation 

As stated in section 3.2, the test functions, v, of the LWF are chosen based on the 
weighted residual (W-R) method being used. Two prominent W-R methods, namely the 
Galerkin and Petrov-Galerkin methods, and their application to beam problems will now 
be discussed. 

In previous literature, a generalized moving least squares (GMLS) interpolation 
scheme was used to develop a Galerkin formulation for solving beam problems (Atluri et 
al . , 1999 ). The trial and test functions in the meshless Galerkin formulation for beam 
problems are chosen to be identical, i.e., Xj = ¥ j ■ This formulation showed 

discontinuities (“scissors” ) at the boundaries of the supports of the trial functions in the 
local sub-domain of the test function (Atluri et al., 1999). Due to these scissors, 
elaborate numerical integration schemes were needed to integrate the weak form 
accurately. The domain of dependence (Q s ) was subdivided into subregions dependent 
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upon where the support domains ended within the Q. s . In each of these subregions, a 10- 
point Gaussian quadrature was used to integrate the weak from accurately. 

In the current work, the Petrov-Galerkin method is used, i.e., the test functions 
are chosen to be distinctly different from the shape functions y/j {% / ^ f// , ) . Recall the 


weight functions chosen as test function components in Eqs. (3.4.6 - 3.4.8), repeated here 
for convenience. 


*# (W) OOH 


1 ~\di/Ro 
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if 0 < dj < R 0 
if dj > R 0 


with di = ||x -Xi\\ and (3= 2,3, and 4, 
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\ R oJ 


+ 2 
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\ R o J 


if 0 < dj < R 0 
if dj >R 0 , 
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x\ w \x)A 


1-6 
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f d L ' 2 

\ R o j 


+ 8 


'i^ 


\ R 0 J 


-3 


'id* 


\ R 0 J 


if 0 < dj < R 0 


if dj > R 0 , 


(3.6.1) 


(3.6.2) 


(3.6.3) 


where R 0 is a user-defined parameter. 

The derivatives of the test functions can be evaluated at the center {dj/R 0 = 0 ) and 
at the end points (di/R 0 = 1) as 


3 W ° r (w) 

3x M( > ' 


d; 


= 0 


\ R 0 J 


fa'" 1 1 


( d- 3 

u i _ j 

K R o J 


= 0 ; Wq > 1 and 

= 0 ; wj > 1 . 


(3.6.4) 
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y 

The test functions are then C continuous up to the order y where y = m i n( m<,. ni\) (see 
Atluri and Shen, 2002). With these definitions, the test functions from Eq. (3.6.1 ) with 
/3 = 2, 3, and 4 are C 1 , C 1 , and C 3 continuous, respectively. Similarly, the spline functions 
from Eqs. (3.6.2 and 3.6.3) are C 1 continuous. As pointed out previously, the lengths R, 
andi? 0 in Eqs. (3.3.19 - 3.3.21 and 3.6.1 -3.6.3) are user-controlled parameters in the 
numerical implementation of the MLPG method. 

To evaluate the validity of the MLPG method and the usefulness of each of the 
trial and test functions derived from Eqs. (3.3.19 - 3.3.21 and 3.6.1 - 3.6.3), the MLPG 
method is applied to various patch tests and mixed boundary value beam problems in 
Chapter 4. 

3.7 Concluding Remarks 

The MLPG formulation was developed for bending of beams - C 1 problems. A 
local weak form (LWF) was developed from the classical weighted-residual form of the 
fourth order governing differential equation. The moving least squares interpolation 
scheme was generalized to include derivatives. These generalized moving least squares 
(GMLS ) approximations were used as the trial functions. The test functions were chosen 
from a different space than the trial functions as combinations of simple weight functions 
and their derivatives. This choice of test functions makes the method a Petrov-Galerkin 
method. Substitution of these trial and test functions into the LWF yielded a system of 
algebraic equations. Stiffness matrices were found to be unsymmetric and banded. Tire 
continuity of the test functions was also discussed. In Chapter 4, several numerical 
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examples are considered to study the effectiveness of the MLPG method for beam 
problems. 
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Chapter 4: Numerical Examples 


Several numerical examples are used to study the effectiveness of the MLPG 
method for beam problems. For the examples presented, a beam of constant flexural 
rigidity El and a length of 4/ is considered. The length 4/ was specifically chosen to 
avoid scaling by unity. Six models with 5, 9, 1 7, 33, 65, and 129 nodes uniformly 
distributed along the length of the beam are considered. Figure 4.0.1 shows a typical 1 7- 
node model. 



Figure 4.0.1: A 17- node model of the beam 


The distances between the nodes (Ax/I) in these models are 1.0, 0.5, 0.25, 0.125, 0.0625, 
and 0.03 125 for the 5-, 9-, 17-, 33-, 65-, and 129-node models, respectively. Three types 
of basis functions, quadratic basis (1, x, x 2 ), cubic basis (1, x, x 2 , x 3 ), and quartic basis (1, 
x, x 2 , x 3 , x 4 ) are used. As mentioned in Chapter 3, the system matrices for the MLPG 
algorithm are of the form (Eq. 3.5.3): 


^ (node) ^ +K (bdry)^_^. (node) _^.(bdry) 


(4.0.1) 


where the superscript “bdry” denotes boundary. These matrices are developed with the 
previously mentioned parameters. Problems studied in this chapter include both patch 
test and mixed boundary value problems. First, simple patch test problems are studied 
wherein a local coordinate approach is developed to improve the accuracy of the method. 
Error norms of the patch tests for both the global and local methods are compared to 


demonstrate the validity of the local approach. Next, general rales of thumb for choosing 


89 



the various user-defined problem parameters are discussed. Then, several mixed 
boundary value problems are worked. Finally, the method is extended to continuous 
beams, and an example problem is studied. As will be demonstrated, the MLPG method 
for beam problems yields excellent results for both primary and secondary variables 
without the need for elaborate post-processing techniques. 

4.1 Patch Tests 

The MLPG formulation for C 1 problems was evaluated by applying the 
formulation to simple patch-test problems. The problems considered were 

1 . w(x) = Cq, 6 = — = 0; Rigid body translation 

dx 

2. w(x) = C\X, 0=c\\ Rigid body rotation (4.1.1) 

2 

3. w(x)=c 2 * 12, 6 = C2X ; Constant -curvature condition 

where Co, C\, and c 2 are arbitrary constants. The third patch test could be looked upon as 
the problem of a cantilever beam with a moment, M=EI(d 2 w/dx 2 )= EJc 2 , applied at x=4/. 
These three problems are depicted in Figure 4.1.1. All three of these problems satisfy the 
governing differential equation exactly and as such represent three simple exact solution 
problems. The deflection w and the slope corresponding to problems 1, 2, and 3 were 
prescribed as essential boundary conditions (EBCs) at x=0 and x=4/. With these EBCs, 
the beam problem was solved using the MLPG method. If the MLPG method recovers 
the exact solution at all the interior nodes and at every arbitrary point of the beam, then 
the MLPG method passes the patch test. 
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• -r- 





(a) Rigid body translation 



(b) Rigid body rotation 



l 4/ l 

(c) Constant-curvature condition 

Figure 4.1.1: Patch tests for beam problems 


In preliminary evaluations, the %\ w \x) test function component in the MLPG 
weak form was chosen as 




i - k 2 /*„ 2 


0 


4 


if 0 < dj < R 0 
if dt>R 0 , 


(4-1.2) 
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where (!■ = || x-x t ||. The weight functions X f (x) used to constnict the trial functions 


were 


Xj(x) = 


l-\dj/Rj 


if 0 <dj<Rj 


if dj>Rj, 


(4.1.3a) 


and 


Xj(x) = 


1-6 


f dj A 


\ R n 


+ 8 


f *j' 


\ R j j 




\ R n 


if 0 < dj < Rj 


if d j > Rj . 


(4.1.3b) 


where d,- = II x - x t II. (Recall from section 3.3 that x = x - x .■ is used in the GMLS 


approximation to construct the trial functions in the local neighborhood x of x. Thus, 
dj = \\x - xj\\ could also be written as dj = ||x|| . } The term (R 0 /I) in the test functions (Eq. 

4.1.2) in each of these six models was different and chosen equal to (2 Ax). The (Rj/l) in 
Eqs. (4.1 .3) were chosen to be (Rj // = 3 .5 ) for the 5-, 9-, and 1 7- node models and (Rj // 
= 16Ax) for the 33-, 65-, and 129- node models. The ranges of (R 0 / 1) and (. Rj / 1) are 
discussed later in this chapter, in section 4.4. 

The displacement vectors {d} that correspond to each of the conditions in Eq. 
(4.1.1) (and in the absence of any other loading) when used in Eq. (4.0.1) should result in 
a null right-hand vector if the K (node) is evaluated exactly. In general, the product results 
in a residual {r (vector as 

K (n ° de) {d} = {r}. (4.1.4) 

Each of the components of the vector { r } is nearly equal to machine zero if K (node) is 
evaluated accurately. To quantify the residual, an error norm of { r } is computed as 
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1 V* 2 

1 \N d ]:' k 


(4.1.5) 


where r k is the k component of the vector {r} in Eq. (4.1.4), and Nd is the degrees of 
freedom in the model. 

Table 4.1.1: Error norm llElli of the residuals for six models and for two basis functions 


Number of 

W= 

-Co 

w= 

C\X 

iv = c 2 x 2 /2 

nodes in 
the model 

Quadratic 

Basis 

Cubic 

Basis 

Quadratic 

Basis 

Cubic 

Basis 

Quadratic 

Basis 

Cubic 

Basis 

5* 

0.5040e-14 

0.1278e-12 

0.2099e-14 

0.4547e-13 

0.5733e-14 

0.9196e-13 

9 * 

0.75 1 5e- 1 3 

0.1496e-l 1 

0.2362e-13 

0.5514e-12 

0.332 le-13 

0.9747e-12 

17* 

0.2774e-10 

0.821 le-10 

0.1 109e-10 

0.3067e-10 

0. 1 582e- 1 0 

0.5352e-10 


0.3609e-9 

0.1691e-6 

0.1796e-4 


0.1062e-5 

0.1435e-2 


0.1266e-9 

0.7735e-7 


0.5599e+0 0.8154e-5 


0.4479e-6 

0.5855e-3 

0.2269e+0 


0.2587e-10 0.9057e-6 


0.1726e-6 
0.1 794e-4 


0.1 193e-2 
0.4154e+0 


* Rj/ 1 = 3.5 


Table 4.1.1 presents the error norm ||E||i for the three conditions in Eq. (4.1.1) 
when the weight function in Eq. (4. 1 ,3b) was used. (Similar results were obtained when 
the weight function in Eq. (4.1.3a) was used. ) As seen from the table, the error norm 
||is||i deteriorates with model refinement and for higher order basis. Closer examination 
of the residuals for each of the six models showed that the residuals were of machine 
accuracy for nodes near the origin while the residuals were largest at nodes farthest from 
the origin. This observation was confirmed by running different cases with the origin at 
different locations along the length of the beam. The origin was moved to the center of 
the beam so that the domain Q became -21 < x <21 . The computed error norms |E||i of 
Table 4.1.1 were best at the center of the beam (at x = 0 ) and inferior at the two endpoints 
x = -21 and x = 21. The origin was then moved to the right end of the beam, i.e., 

- 4/ < x < 0 . The error norms |Ej|i were found to be best at the right end of the beam 
(x = 0 ) and inferior at the left end of the beam (x = -47). In fact, the same error norms of 


93 
















0.1794e-4 and 0.4154e+0 for the 129-node model and w = c 2 x 2 / 2 were observed at 
(x = -41) when the origin was moved to the right end of the beam. Also, the residuals 
were largest for the models with the largest number of nodes. This behavior is counter- 
intuitive. 

Closer scrutiny of computations showed that the numerical values of the shape 
functions for nodes that are systematically located about the center of the beam ( for 
example, nodes 3 and 15, 2 and 16, and 1 and 17 in the 17-node model of Figure 4.0.1) 
are not exactly identical as expected. These differences increased with model refinement 
and when a higher order basis was used. The reason for this behavior is explained in 
section 4.2. 

4.2 Local Coordinate Approach 

In the MLS interpolation, the basis functions are in terms of the global coordinate 
x. The [A] matrix thus formed using this basis is generally of the form (see Atluri et al. , 
1999, Eq. 16) 

[A] = ^ \l k (x) p p T + A k (x) p x ■ pi } (4.2.1) 

k = 1 

where x = x-Xj, and Mare the number of nodes in the domain of definition of node j 

for which the [A] matrix is being computed. (For convenience in presentation, the [A] 
matrices thus formed will be referred to as the global method. ) As the order of the 
polynomial basis increases, the conditioning of the [A] matrix deteriorates. For example, 
the matrix [A] will have terms like 1 , x 2 , x 4 , x 6 on the diagonal for a cubic basis function. 
The [A] matrices for nodes near the origin and the [A] matrices for nodes farthest from 
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the origin will be different. The conditioning is worse for [A] matrices for nodes farthest 
from the origin. This explains the differences in the error norms observed in Table 4.1.1. 
The error norms in Table 4.1.1 can be improved by using higher precision computations 
or inversion routines. However, a much simpler alternative to improve the accuracy is 
discussed. 

The conditioning of the [A] matrix can be considerably improved if the MLS 
approximation is defined not in terms of a global basis, but rather in terms of a local 
basis. Figure 4.2.1 shows two identical shape functions, one centered at node j, and the 
other centered at node e. 



2 Rj 

Figure 4.2.1: Local coordinate definitions 

The global approximation for 


2 Rj 


w(x) = p (x)a(x) 


2 m — 1 

■ci\+a 2 X + <? 3 X +... + a m x 


(4.2.2) 


can be rewritten in the neighborhood of node j , recognizing that x = Xj + E, where ^ is a 
local coordinate measured from node j, as 
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w(x) = ti\ + a 2 [xj + ^)+ «3 {xj +£;) 2 + . . . 

= + ci 2 x j + ai,Xj + . . .)+ ( 0 2 + 2 a^Xj + • • + (#3 + . . . ) t/' 


(4.2.3) 


-b\+ b 2 % + b 3 ^ + . . . + 


where bt,i = 1,2,..., m are the new undetermined coefficients in the MLS 


approximation. (A similar local coordinate transformation can be affected for node e in 
Figure 4.2.1 as x = x e + <g.) The [A] matrix then is computed in a similar manner as in Eq. 


(4.2.1), but with 


pV)= 1, l Z 2 , ... 1 


and pI(£> = fo, 1 , 2 £ 3^ 2 , ... 


(m-2) 


(4.2.4) 


as A( )*( ) . 
ax dg 

The local coordinate approach was implemented in the evaluation of the shape 
functions and their derivatives for all the nodes in the six MLPG models of the beam. 


Table 4.2.1 compares the condition numbers of the [A] matrices at various locations on 
the beam using global and local coordinate methods. The condition numbers were 
evaluated using routines available in NAPACK and the procedure outlined in Bathe, 

1996 and Chapra and Canale, 1988. A brief review of condition numbers is presented in 
Appendix B. When the global coordinate method was used, the condition numbers of the 
[A] matrices for nodes farthest from the origin were much larger (suggesting poor 
conditioning ) than the nodes closest to the origin. The conditioning numbers of the [A] 
matrices vastly improved when the local coordinate method was used, clearly 
demonstrating the advantages of the local coordinate method. 
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Table 4.2.1: Comparison of the condition numbers of the [A] matrices at various 
locations on the beam using global and local coordinate methods 


Location 
on the 
beam (x/4 1) 

Number of nodes in the model 

5* 

9 * 

17* 

33 

65 

129 

Global Method Conditioning Number 

0.0 

0.631e+3 

0.106e+4 

0.930e+3 

0.271 e+3 

0.267e+3 

0.189e+4 


0.23 le+5 

0.268e+5 

0.272e+5 

0.785e+5 

0.904e+6 

0.13 le +8 


0.914e+6 

0.771e+6 

0.127e+7 

0.422e+8 

0.153e+10 

0.365e+l 1 


Local 

1 Method Conditioning Number 

0.0 

0.634e+3 

0.106e+4 



0.267e+3 

0.189e+4 


0.478e+3 

0.496e+2 

0.41 le+2 

0.11 le +2 

0.153e+2 

0.141 e+3 


0.632e+3 

0.106e+4 

0.930e+3 

0.271 e+3 

0.267e+3 

0.189e+4 


* Rj/l= 3.5 


The error norms shown in Table 4.1.1 were recomputed and the results are 


presented in Table 4.2.2. As expected, all models and the quadratic and cubic basis 
functions produced error norms close to machine accuracy, suggesting that the local 
coordinate approach produces a significant increase in accuracy compared to the global 


coordinate approach. 


Table 4.2.2: Error norm H^Hi of the residuals computed with the 


Number of 
nodes in 
the model 

w=Co 

W=C\X 

W=C 2^12 

Quadratic 

Basis 

Cubic 

Basis 

Quadratic 

Basis 

Cubic 

Basis 

Quadratic 

Basis 

Cubic 

Basis 

5* 

0. 1 1 73e- 14 

0.3500e-13 

0.2342e-15 

0. 1 20 1 e- 1 4 

0.3174e-14 

0.3853e-13 

9 * 

0.2521e-13 

0.4900e-13 

0.8357e-14 

0. 1 699e-l 3 

0.3659e-13 

0.4146e-13 

17* 

0.1392e-12 

0.2169e-12 

0.4764e-13 

0.1680e-12 

0.2126e-12 

0.8124e-12 

33 

0.4389e-12 

0. 1 390e- 1 1 

0.1876e-12 

0.5060e-12 

0.4084e-12 

0.2 1 83e-l 1 

65 

0.4 1 96e- 1 1 

0.3890e-l 1 

0.1 142e-l 1 

0. 1 879e- 1 1 

0.2548e-l 1 

0.5930e-l 1 

129 

0.4029e-10 

0.2778e-10 

0.1240e-10 

0.8 1 91 e -1 1 

0.2400e-10 

0.2166e-10 


ocal coordinate approach 


■Ri/l = 3.5 


In the global MLPG implementation, the [A] matrix is calculated and inverted at 


every node in the model. When using the local coordinate methodology with uniform 


nodal spacing, the shape functions are exactly identical for nodes whose R f places the 
entire shape function in the interior of the domain of the problem. Hence, for those nodes 
the [A] matrices are identical. As such, considerable reduction in computational effort 
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and cost may be achieved by the proposed local coordinate approach. The local 
coordinate approach is used to evaluate the shape functions and their derivatives for all 
numerical examples in the remainder of this report. 

4.3 Patch Tests Revisited 

The three patch test problems, 

1 . w(x) = Cq, 9 = — = 0; Rigid body translation 

dx 

2. w(x) = qx, 9=c\\ Rigid body rotation (4.3.1) 

9 

3. w(x) = C 2 X 12, 9-C2X\ Constant -curvature condition 

outlined in section 4.1 are now studied. 

For a displacement of w(x) = c 0 and C\X units, the rigid body conditions were 
modeled with boundary conditions 

Translation: Rotation: (4.3.2) 

w|x = o=c 0 w/ x =41=Co w/ x = 0=0 w/ x = 4l - Ac x l 

9 |a = o=0 9 |.v = 4 / = 0 9\ x = q=C\ 9\ x = 4 i=C\ 

Since the exact solutions are constant and linear in x, respectively, the MLPG method 

developed with a quadratic or higher order basis function must reproduce the solutions 

exactly. (A linear basis cannot be used as the LWF requires second derivatives of the 

trial functions.) As expected, the algorithm reproduces the exact solutions for w> and 6? to 

machine accuracy for both rigid body modes at the nodes and at any arbitrary point in the 

beam. 

For the constant - curvature condition, w = c 2 x 2 /2, the problem was modeled with 

EBCs 
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(4.3.3) 


w/x = 0 - 0 w/x = 41 - 8C2 / 2 

0\x = o= 0 d\ x = 4i = 4 c 2 1 

Since the exact solution is quadratic in x, the MLPG method developed with a quadratic 
or higher order basis function must reproduce the solution exactly. As expected, the 
algorithm reproduced the exact solution for the primary variables to machine accuracy at 
all nodes and at any arbitrary point in the beam. 

The above analyses were repeated with each of the test function components in 
Eqs. (3.4.6 - 3.4.8). The MLPG method reproduced exact solutions to machine accuracy, 
thus passing all the patch tests. 

4.4 Problem Parameters 

As mentioned previously, the parameters ( R 0 / 1 ) and (Rj / 1) in the MLPG method 
are user-controlled. Ranges of values of these parameters were studied, and a general 
rule of thumb was established. The previously mentioned lengths (R 0 / 1 = 2 Ax) and 
(Rj/ 1 = 8 Ax) were used at all nodes of an A-node model, except at node 2 and node A- 1 
(see Figure 4.4. 1). For these nodes, (R 0 / 1 = Ax), was used to ensure a symmetric Q. s and 
account for the terms of Eqs. (3.2.18) where, with (R 0 /l = 2Ax) for nodes 2 and A-l, 

0 < V 2 < 1 and 0 < Vn- i < 1 . Note that with these assignments of (R 0 / 1) the test functions 
for all interior nodes have symmetric Q. s configurations. As shown in Figure 4.4.1, no 
asymmetry is introduced at nodes 1 and N as exactly half of their test functions are used. 
When these symmetries are violated, the MLPG method requires additional terms as 
discussed in section 3.2. When these terms are accounted for, the MLPG method passes 
the patch tests. 
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Figure 4.4.1: Local sub-domain, Q s , definitions for various nodes 

As the models are refined, the value of (At //) decreases and thus the size of Q s 
and the extent of the trial functions also decrease. For finer models, i.e. for the 33-, 65-, 
and 129-node models, when 8 Ax' <Rj 11 < 16Ax the MLPG method yielded very 

accurate results. However, when Rj/ 1 > 16 Ax' the MLPG method failed the patch tests 
for these models. Figure 4.4.2 shows the results of the rigid body rotation problem for 
these two cases. When R,-/l > 16Ax the trial function is too diffused and the size of Q. s 
(Ro / 1 = 2 Ax) is too small in comparison to (Rj / 1). The combination of small Q s size and 
large (Rj / 1) are apparently incompatible. While the finer models performed well over a 
large range of (Rj / 1), the coarser models performed well in a much smaller range of (Rj / 
1). For good performance, (Rj / 1) needed to be approximately 8Ax but less than 98% of 
the total beam length. 
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MLPG 
Exact 


(a) MLPG and exact solutions when 8Ax < Rj/l < 16Ax 



— © — MLPG, satisfactory 
value of Rj 

— a— mlpg, unsatisfactory 
value of Rj 


(b) MLPG solutions when R -/1 > 16Ax 


Figure 4.4.2: Rigid body rotation - Comparison of results for 
different extents of trial functions 


4.5 Mixed Boundary Value Problems 

The MLPG method was applied to beam problems with mixed boundary 
conditions. Recall from section 3.3 that the prescription of displacement and shear and 
slope and moment are mutually disjoint. When displacement is prescribed, the shear 
cannot be prescribed. Likewise, when slope is prescribed, the moment cannot be 
prescribed. 
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4.5.1 Cantilever beam with concentrated moment at the free end 

The first problem considered was a cantilever beam with a concentrated moment 
at the free end (i.e. M - Mo atx = 41, see Figure 4.5.1). 



Figure 4.5.1: Cantilever beam with concentrated 
moment at the free end 

The exact solution for this problem is w = Mox 2 / 2EI and 0 = Mox / El. For all trial 
functions considered, the MLPG algorithm reproduced the exact solution when the test 
function components in Eq. (3.4.6) with J3= 3 and 4 and when the 4-term spline function 
in Eq. (3.4.8) were used. In contrast, the algorithm failed to reproduce the exact solution 
when the test function component in Eq. (3.4.6) with J3= 2 and the 3-term spline function 

of Eq. (3.4.7) were used. This example suggests that ) test function components with 

at least C 1 continuity and with m x > 2 (see Eq. 3.6.4) are required for the MLPG 
algorithm for beam problems. 

4.5.2 Cantilever beam with tip load 

The second problem considered was a cantilever beam with a tip load (See Figure 
4.5.2). Since the exact solution for this problem is cubic in terms of the x-coordinate of 
the beam, all six models with a cubic basis function and a test function with C 1 continuity 
and with m x > 2 reproduced the exact solution to machine accuracy. 
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4.5.3 Simply supported beam subjected to uniformly distributed load 

The third problem considered was a simply supported beam subjected to a 


uniformly distributed load (see Figure 4.5.3). 





Figure 4.5.3: Simply supported beam subjected to 
a uniformly distributed load 


The exact solution for this problem is given by 


q ( ~ T 3 4 r 3 dw q ( . T 2 . 3 r 3 

w = — - — -2 Lx + x +Lx\ , — = — - — -6 Lx +4x +L 
24 El \ ) dx 24EI V 


where L = 41. Using symmetry, half of the beam was modeled. Since the exact solution 
for this problem is quartic in terms of the x-coordinate of the beam, the MLPG method 


with a cubic basis function did not reproduce the exact solution. Error norms defined as 
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were computed at g uniformly spaced points along the beam. A value of g = 200 was 
used. The norms \\E W \\ 2 and \\E M \\ 2 are presented in Table 4.5.1 . As expected, all models 
(N>9) yielded accurate solutions (within 4% for w and M). As the number of nodes in 
the models were increased, the ||is w ||2 norm changed marginally, suggesting the same 
accuracy in the solutions for the various models. Also, the \\E M \\ 2 norm was of the same 
order as the ||£’ w ||2 norm, suggesting the same accuracy for the primary and the secondary 
variables. To obtain acceptable results using a Galerkin formulation ( Atluri et al., 1999), 
El s would have to be subdivided into sub-regions within which, for example, a 10-point 
Gaussian quadrature would be used to perform the integrations (see section 3.6). The 
numbers in Table 4.5.1 were computed via a Petrov-Galerkin formulation using a 20- 
point Gaussian integration in each of the single compact support domains Q. s . When the 
order of the basis function was increased to quartic, the MLPG method reproduced the 
exact solutions (for w, 6 \ M, and V) to machine accuracy. 

Table 4,5.1: Error norm ||E|| 2 for a simply supported beam subjected to a uniformly 
distributed load with cubic basis used in the MLPG method. (Trial function using 
Eq. (3.3.19) with oc=3 and test function using Eq.(3.4.6) with J3=4.) 


Error norm 

Number of nodes in the moc 

lei 

5* 

9* 

17 t 

33 t 

65 t 

129 t 

E w 

2 

0.1662e-l 

0.1306e-2 

0.4573e-2 

0.3829e-l 

0.1742e-l 

0.2368e-l 

Em 

2 

0.2774e+0 

0. 1 057e- 1 

0.1704e-l 

0.3680e-l 

0.1763e-l 

0.2340e-l 


* Rj / 1 = 3.5,^ Rj / 1= 8Ax 


The previously discussed problem demonstrates an interesting phenomenon. 

When the order of the basis function equals the order of the exact solution, the previously 
discussed 8-point Gaussian quadrature in a single Ei s is sufficient to integrate the weak 
form very accurately. However, when the order of the basis function is less than the 
order of the exact solution, a higher order integration rule (such as a 20-point Gaussian 
integration) is needed to obtain accurate results. For problems with complicated loading, 
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where exact solutions are not known, the order of the basis function can easily be 
increased until convergence of the solution is achieved. 

The problem of the simply supported beam subjected to a uniformly distributed 
load was modeled next using the full beam with non-uniform nodal spacing shown in 
Figure 4.5.4. 



(— > x 

Figure 4.5.4: A 19-node model with unequally spaced nodes 

This model was generated by randomly placing nodes in the region 0 < x < 21 and 
symmetrically replicating these nodes in the region 21 <x <41 . The order of the basis 
function was increased to quartic. The MLPG and exact solutions for deflection, 
moment, and shear are presented in Figure 4.5.5. 



0 0.5 1 


x 1 41 

(a) Deflection 

Figure 4.5.5: MLPG and exact solutions for a simply supported beam 
subjected to a uniformly distributed load 
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x / 4/ 

(b) Moment 

1 

0.5 

V(x) 

^max 0 
-0.5 



Figure 4.5.5 Concluded: MLPG and exact solutions for a simply 
supported beam subjected to a uniformly distributed load 

As expected, the MLPG method reproduced the exact solutions to machine accuracy for 
both the primary and secondary variables despite the nodal arrangement. 

4.5.4 Simply supported beam subjected to a central concentrated load 

The fourth problem considered was a simply supported beam subjected to a 
central concentrated load (see Figure 4.5.6). 
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p 



Figure 4.5.6: Simply supported beam subjected to a 
central concentrated load 
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(Rj / 1 ) was chosen as (R, 1 1) = 8Ax for the 33-, 65-, and 129-node models, and as (Rj / 1) = 
3.5 for the 5-, 9-, and 17-node models. For the full representation of the beam, the values 
of (Rj / 1) are noted in Table 4.5.2. The exact solution for the deflection under the load is 
given by 


PIS' _P(41 ) 3 
48 EI~ 48 El ‘ 


The exact solutions for the slopes at the end points are given by 



PL _ 
16 El 


pm 2 

\6EI 


and 0 \ x=4 /=- 


PL 

\6EI 


P(W 

\6EI 


(4.5.5) 


(4.5.6) 


For the symmetric representation of the beam, the boundary conditions shown in 
Figure 4.5.7 were used. As expected, the MLPG method reproduced the exact solutions 
for all models at all nodes and at all interior points of the beam. 


w = 0 



V = -PI2 


M = 0 9=0 


Figure 4.5.7: Symmetric representation of a simply 
supported beam subjected to a central concentrated load 


The concentrated load at the center of the beam is expected to cause difficulty 
when the full beam is considered. As such the full beam is modeled to study the 
performance of the MLPG method. External loads contribute to the f tnode) of Eq. (4.0.1) 
(see Eq. 3.5.4f), repeated here: 
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j.(node) 



(4.5.7) 


In numerical implementation, if a concentrated load P is applied at node p, the integrals 
of Eq. (4.5.7) are evaluated with the dirac delta function as 


„(node) 


\zi W> fdx 


\x\ W) f 3(x = x p ) dx 



A 


q(0 


!T’’V>v 

Kv * 

> — < 

\x\ 6 > f Six = x p )dx 

> — < 

xf'lXpiP 

Up J 


Q (0 

l J 




(4.5.8) 


To evaluate the f (node) , all the nodes in the domain of influence of node p need to be 
examined. The value of each test function, v„ in the domain of influence of node p is 
evaluated at node p. As the values of the test functions, v t , are computed as ( see Eq. 


3.4.3) 


v/ (x p ) = /4 W V ; (W) (x p ) + p { i d) x\ d) (x p ) , (4.5 .9) 

the corresponding and xf ] contribute to the f (node) as in Eq. (4.5.8). 

For the full representation of the beam, the MLPG values of ww and 9 ma for 
each of the models studied are presented in Table 4.5.2. 


Table 4.5.2: MLPG values of deflection and slope for models with various nodal arrangements 



Numl 

3er of nodes in the model; (/?, / 1 ) 

5; (3 Ax) 

9; (4 Ax) 

17; (6 Ax) 

33; (8 Ax) 

65; (8 Ax) 

129; (8Ax) 

^(max) / W(max) exact 

0.9746 

1.0882 

1.0368 

0.9982 

0.9992 

1.0120 

@ (max)/ ^(max) exact 

0.3717 

1.1003 

1.0380 

1.0012 

0.9975 

1.0126 


The MLPG and exact solutions for deflection and moment of the 65 -node model are 
compared in Figures 4.5.8. 
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O MLPG 
— exact 


(a) Deflection 


l 


M(x) 

^max 


0 

0 0.5 1 

x/AJ 

(b) Moment 

Figure 4.5.8: MLPG and exact solutions for a simply supported beam 
with a central concentrated load 

These figures and the results presented in Table 4.5.2 demonstrate that the MLPG method 
yields excellent results for both primary and secondary variables. These results were 
obtained without the use of elaborate post-processing techniques. As the number of 
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nodes was increased from 17 to 129, the accuracy of the solutions did not appreciably 
change, suggesting that a 17-node or 33-node model is sufficient to obtain an accurate 
solution. The MLPG method apparently handled the discontinuity caused by the central 
concentrated load well. 

4.6 Continuous Beams 

The MLPG method was then applied to a continuous beam problem to evaluate its 
effectiveness. A continuous beam with one additional support along the interior of the 
beam (shown in Figure 4.6.1) is considered. 

z 

( 1 

1 [ i [ ; r 

>h 

l l 

Figure 4.6.1: Continuous beam subjected to a uniformly distributed load 

In applying the MLPG method to continuous beams, an additional penalty term, 

a c [{w-w)v \ rf , (4.6.1) 

where a c is the penalty parameter to enforce the continuous beam boundary condition, is 
added to the weak form. The weak form of the governing differential equation then 
becomes 




-> x 
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) = EI [ d Wd V dx- [f vdx + aJiw-w)v\ T + a e 

J dx 2 dx 2 J w LW* )dx\ Te 


3 i r 2 

w 1 „,(/ w dv nr d w 

+ a c [{ W-M'jvjr +«X vEI T- -»x ~r EI r 

' A- dx dx 2 


(4.6.2) 


As in section 3.2, utilizing the boundary condition subsets, 


r,nr w , r s nr 0 , 
r v n ry , and r, n r M . 


(4.6.3) 


leads to the local weak form (LWF) for continuous beam problems: 


0 = El \ d l dx - f / v dx 

J dx dx J 


+ a w [(w - ^)v]r w +a e + a c - ^>]r sc 

r r 3 n r 2 

[\7 1 jU dv rr, d w r . d wdv 

-«xlvvk„ ~n x M— +n x El — v -n x EI—~— 


(4.6.4) 


In comparing Eq. (4.6.4) with the LWF developed in section 3.2 (see Eq. 3.2.15), it is 
noted that the only difference in the two equations is the teim a c [(w- M^)v] r . 

Therefore, the LWF of Eq. (4.6.4) can be used for all beam problems worked in this 
report; when no continuous beam boundary conditions are present, a c = 0, and Eq. (4.6.4) 
becomes Eq. (3.2.15). Following the development of section 3.5, the MLPG equations 


K (node)^ + K <bdry)^ _ f (node) _ f (bdry) _ 


(4.6.5) 


where K (node) and f" odc ' remain as Eqs. (3.5.4b, 3.5.4d, and 3.5.4f), and the expressions 
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(a) Deflection and Slope 


Figure 4.6.2: MLPG and exact solutions for primary and secondary 
variables of a continuous beam subjected to a uniformly distributed load 
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Exact 

MLPG 


(b) Moment 

Figure 4.6.2 Concluded: MLPG and exact solutions for primary and secondary 
variables of a continuous beam subjected to a uniformly distributed load 
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Chapter 5: Concluding Remarks 


The Meshless Local Petrov-Galerkin (MLPG) method has been implemented for 
2-D potential and elasticity problems. In this report, the method was implemented and 
studied for 1-D C° problems and further developed for bending of beams - C 1 problems. 
The following conclusions are drawn from the work presented in this report: 

• The MLPG method yields accurate solutions for C° and C 1 problems. 

• The MLPG method yields continuous secondary variables as demonstrated by the bar 
and beam problems studied. 

• A local coordinate approach is developed and validated for improving the 
conditioning of the [A] matrix that is needed to evaluate the trial functions for the 
beam problems. 

• For beam problems, the Petrov-Galerkin approach is preferable over the Galerkin 
approach. 

• Reasonable ranges of several parameters are required to obtain good results. The 
ranges of these parameters suggest the robustness of this method. 

Each of these conclusions is discussed below. 

5.1 Accurate Solutions by the MLPG Method 

As discussed in Chapter 1, for any new method to compete with the Finite 
Element Method, the new method must retain the advantages of the FEM. This includes, 
most importantly, the ability of the method to yield accurate ( to machine accuracy) 
solutions. (As stated in Chapter 1, “machine accuracy” means that the difference 
between the exact and numerical solutions is of the order of 10" 14 when double precision 
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arithmetic is used.) For all the C° and C 1 problems presented in this report, the MLPG 
method yielded accurate solutions for both the primary and secondary variables. 

5.2 Continuous Secondary Variables 

As discussed in Chapter 1, one of the disadvantages of the FEM is the 
discontinuity of the secondary variables across inter-element boundaries. The 
discontinuities in the secondary variables arise because of the piecewise linear shape 
functions that are used to construct the trial functions. Elaborate post-processing 
techniques are needed to obtain smooth distributions of these secondary variables. In the 
MLPG method, elements are eliminated, and nodes are utilized in the domain of the 
problem. A diffused (i.e., not piecewise linear) trial function such as a moving least 
squares (MLS) interpolation is used. These diffused trial functions are smooth, and 
hence, smooth distributions of the secondary variables are obtained, thus eliminating the 
disadvantage of the FEM. These results were confirmed in Chapter 2 by application of 
the method to a C H problem and in Chapter 3 for C 1 problems. 

5.3 Local Coordinate Approach 

As discussed in Chapters 2 and 3, the trial functions used to approximate the 
solution are formed from shape functions that are developed by a MLS interpolation. 

The formation of these shape functions involves the evaluation of the [A] matrix. The 
[A] matrix is evaluated using the weight functions and the basis functions. The 
conditioning of the [A] matrix is determined by the order of the basis function used. As 
the order of the basis function is increased, the conditioning of the [A] matrix becomes 
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poor, especially for nodes far from the origin, resulting in an inaccurate computation of 
the inverse of [A] that leads to poor quality solutions. To improve the conditioning of the 
[A] matrix, the MLS approximation is defined in terms of a local basis rather than a 
global basis. A comparison of the results of the global and local approaches applied to 
patch test problems as presented in section 4.2 clearly demonstrates that the local 
coordinate approach produces very accurate results in comparison to the global 
coordinate approach. 

5.4 The Petrov-Galerkin Approach 

In the MLPG method, because the trial and test functions are chosen from 
different spaces, the resulting system stiffness matrices are unsymmetric. This could be 
perceived as a disadvantage. However, closer examination of the method reveals that this 
is not a disadvantage. As discussed in section 3.6, when a Galerkin approximation is 
used, the system matrix is symmetric. However, the Galerkin approach results in 
discontinuities that arise at the boundaries of the supports of the trial functions in the 
local sub-domain of the test function. Because of these discontinuities, elaborate 
numerical integration schemes are needed to integrate the weak form accurately. The 
local sub-domain, Q s , of the test function is divided into sub-regions. The endpoints of 
these sub-regions are determined by the ends of the support domains of the trial functions 
that intersect Q s . A 10-point Gaussian quadrature is used in each of the sub-regions to 
accurately integrate the weak form. The sub-region procedure results in large computing 
effort to integrate the weak form in each sub-domain L2 y . Alternately, if a Petrov- 
Galerkin approximation is used, the sub-region integration is not needed. A single higher 
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order integration rule (for example, a 20-point Gaussian) in a single sub-domain is 
sufficient to integrate the weak form accurately. This result was confirmed numerically 
for several examples in Chapter 4. Thus, while the unsymmetry of the stiffness matrices 
in the Petrov-Galerkin method may be construed as a disadvantage, it is far outweighed 
by the computational time and effort saved by the weaker requirements for numerical 
integration. 

5.5 Problem Parameters 

In applications of the MLPG method, several parameters are user-defined. Over 
certain ranges of these parameters, good performance is obtained. The minimum order of 
Gaussian integration required depends on the basis functions and weight functions used. 
Also, extremely high orders of Gaussian integration are unreasonable and unnecessary. 
For the problems worked in this report, numerical experimentation showed that a 20- 
point Gaussian, while not necessary for all simpler problems, was found to integrate the 
weak form accurately. The algorithm performs best when the extents of the test functions 
are in the range Ax' < (R 0 //) < 2 At, where Ac is the nodal spacing between nodes for a 
uniformly distributed nodal arrangement. Similarly, the extents of the trial functions are 
best chosen as 8 A' < (Rj //) < 16 A, but no larger than 98% of the domain of the 

problem (for 1-D problems). 

5.6 Contributions of this Research 

Meshless methods are becoming increasingly popular as evidenced by the large 
amounts of literature on the subject published in the last five years. However, much of 
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the research that is being conducted on meshless methods is on C° problems. In this 
report, one particular meshless method, the Meshless Local Petrov-Galerkin (MLPG) 
method, was extensively studied for C 1 problems. At the time this research was 
conducted, the literature available on the MLPG method for beam problems utilized a 
Galerkin approach. In this report, a Petrov-Galerkin approach was implemented and 
shown to be far superior to the previously available Galerkin approach (see section 3.6). 
Additionally, four major contributions of this report to the general field of meshless 
methods are: (1) a local coordinate approach to be used in the formulation of the trial 
functions was proposed and validated in section 4.2, (2) the performance of several test 
functions (presented in section 3.4) was studied, and well-defined continuity 
requirements for prospective test functions were determined, (3) application of the 
method to a problem with load discontinuity was demonstrated in section 4.5, and (4) 
application of the method to continuous beams was also demonstrated (section 4.6). The 
following publications came out as a result of the work performed in this report: 

I. I. S. Raju and D. R. Phillips (2002): “A Local Coordinate Approach in the MLPG 
Method for Beam Problems,” NASA TM-2002-21 1463, and 

II. I. S. Raju and D. R. Phillips (2002): “A Meshless Local Petrov-Galerkin Method for 
Euler-Bemoulli Beam Problems,” Proceedings of the ICES ’02 conference, Reno, 
Nevada, July 3 1 - August 2, 2002, Paper No. 1 39. 
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5.7 Suggestions for Future Work 

The MLPG method is still in the early stages of its development. More work 

needs to be performed before the method can reasonably compete with the FEM. The 

following are the next steps to extend the research conducted in this report. 

• The method could be extended to Timoshenko beam problems. (First order shear 
deformation is accounted for. The assumption of normals before deformation 
remaining normal after deformation is relaxed to “normals remain straight but need 
not be normal after deformation.”) 

• The method could be extended to two dimensions for plate bending. 

• The method needs to be modified and applied to study built-up structures. 

• The method could be studied for shell analysis. 

• The method relies very heavily on the user-defined parameters, namely R 0 , and Rj. 
More research should be done to determine more robust ranges of these parameters so 
the method consistently obtains good results. 
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Appendix A: Computation of Derivatives of Shape Functions 

This appendix presents a detailed derivation of the derivatives of the shape 
functions used in the MLPG method. Section 1 presents the derivatives of the shape 
functions for C° problems. Section 2 discusses the [B] matrix for C 1 problems. Section 3 
presents the derivatives of the shape functions for C 1 problems. 


A.l C° Problems 

In this section, the derivatives of the shape functions for C° problems are derived 
first in terms of the general spatial coordinates, %*, and then reduced to one dimension. In 
C° problems, the approximations for the solution, u h (x), can be written as 

u h (x ) = p 1 (x)a(x) , (A. 1.1a) 

where p is an rn h order basis functions and a is a vector of undetermined coefficients, and 


’(x) = J^jWuj , 
7=1 


(A. 1.1b) 


where (f>j{x) are shape functions, and u ,• are fictitious nodal values. In Eqs. (A.l .1), x 


represents the spatial coordinates, 
x = [*i *3 1'' • 


(A. 1.2) 


For the local weak form, the derivative of u(x) is needed. First consider the statement of 
u ,l {x) in Eq. (A. 1.1a). Differentiating, 


du 7 (x) T, , r)a(x) dp T (x) 

= P (x)— +-^ a(x). 


dx k 


dx k dx k 


(A. 1.3) 
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In Eq. (A. 1.3), the second term is easy to evaluate; however, evaluation of da (x)/dx k in 
the first term is not straightforward. The evaluation of this derivative requires closer 
examination and is discussed below. 

Consider the equation 


[A] {a} = [B] {h} , 

1) (hi, «)(«,!) 


(A. 1.4) 


in which {u} are fictitious nodal values, and [A] and [B] are easily evaluated from 
weight functions and basis functions using Eqs. (2.2.17 and 2.2.18). Eq. (A. 1.4) can be 
differentiated as 


[A] M + |Ai |a)=[B] Ml + M {ai . 

uXfc OXfc OXfc OXfc 


(A. 1.5) 


In Eq. ( A. 1.5), the fictitious nodal values {u} are not functions of x, and thus the term 
[B] ■ (3 {u}/3x£ ) vanishes. Rearranging Eq. (A. 1.5) one obtains 


[A] M = M {a! _Ml |a| 

OXfc OXfc OXfc 


(A. 1.6) 


This leads to 


9{a 

dxj 


-1 3[B] 


1 =[A] *^{u}-[A] 


dxi. 


■ ia[A] ‘a' 

dx k 


(A. 1.7) 


The vector {a} can be evaluated using Eq. (A. 1.4), 

{a} = [A] -1 [B] {u}. 

(hi, 1) (hi, hi) (m,n)(n, 1) 

Substituting Eq. (A. 1.8) into Eq. (A. 1.7), 

M = [Ar , M, a) _[A]-'|Al [Ar > [ B ]ii ). 

()x/ c dx k dx k 

So, substituting Eqs. (A. 1 .9 and A. 1 .8) into Eq. (A. 1 .3), 


(A. 1.8) 


(A. 1.9) 
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(A. 1.10) 


A^=p t (x, 

OX 7. 


[Ar'ffliinAf'ffliAr'iBnii 


+ 


9p T (x) 


dx t 


dxi- 


[A] \B]{u} 


dx t . 


Now consider the statement of u‘( x) in Eq. (A. 1 . lb). Differentiating, 


du (x) 


Z UU j 

U/x)— + 

7= l v 


du ; dtp Ax) ^ \ dtp Ax) 


dXj ( \ T 7 ' 7 dx k dxj. 1 


=n- 

7=1 


dXl; 


Equating the two expressions for the derivative of u (x), i.e., 
du h (x) I 


dxi 


du h ( x ) 


Eq. A. 1.10 


dXi 


Eq. A. 1.1 1 


leads to 


P 0){u} 


[ A ] 'M-[ A ] '^1[A] '[B] 


dx k 


dx k 

3p 1 (x) r , ,-i 


dxi. 


[A] 1 [B]{u} = ^] 


b d<P,(x) 


7=1 


dxi- 


or. 


<¥_/(*) T rA1 -l rD1 T 
— - — = P ,k [A] [B] + p 

o^k (l,w) (m,m) ( m,n ) (l,w) 


[A] ! [B] ^.-[A] 1 [A] ^ [A] 1 [B] 

(m,m) (Wj „) (m,m) ( m>m ) (m,m) (m,n) 


where 


k = 1,2,3. 

’ dx k 


Finally, 


(A. 1.11) 


(A. 1.12) 


(A. 1.13) 


(A. 1.1 4) 


(A. 1.15) 
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0j,k =YPg,k{\- A \ \v\) +Pg\[M ^A^tA] \b]\. (A. 1.16) 

*=1 'gj J 

Note that in this report, there is only one spatial coordinate, x=x, as 1-D problems are 
considered. As a result, Xk = x, and the partial derivatives become full derivatives: 

TTl 

^A=Z^^( [Arl[B 4 + ^l A r 1 [ B lx-[Ar 1 [A] ;X [A]- 1 [B^ (A. 1.17) 

g=l 

where 

o^Al. 

A.2 |B] Matrix for C 1 Problems 

In section A.l, the [B] matrix is an (m,n) matrix given by Eq. (2.2.18) as 
[B]= P l l . (A.2.1) 

In C 1 problems, the [B] matrix is an (m,2n) matrix given by Eq. (3.3.14) as 

[B]= P l l, pJX, . (A.2. 2) 

Consider the equation 

[A] {a} = [B] {s} . (A.2. 3) 

1 ) (m,2n) (2n,X) 

in which, as in Eq. (3.3.8), {s} is a vector containing {w} and {t } , the fictitious nodal 
values of deflection and slope, respectively. Also, [A] and [B] are easily evaluated from 
weight functions, basis functions, and derivatives of basis functions using Eqs. (3.3.13 
and 3.3.14). The approximations to the solution, w h (x), can be written as in Eq. (A.1.2a) 
as 
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(A.2.4) 


w ,l (x) = p T (x)a(x) . 

Solving for {a} in Eq. (A.2.3) and substituting into Eq. ( A.2.4) yields Eq. (3.3.16), 

m/'(x) = P T (x)[a]- 1 [B] IT], (A.2.5) 

(1,/m) (m,m)(m,2n) l*J 

(2n,l) 


or 


w 


! (x) = p T (x)[A] 1 [B w ]{w}+p T (x)[A] 1 [B ? ] {t} , 

(1,/m) (m,m)(m,n)(n,\) (1,/m) (/m,/k) (m,n)(n,\) 


where 


[B w B,] = 


P l l pjx, 


From Eq. (A.2.6), the shape functions are as Eqs. (3.3.18): 


¥ 


m 

;. vv) (x) = ^p g (x 7 -)[A- 1 P 1 \ 


g=l 

m 


¥ { f } (x) = 5>,(x y H A lp v ^ 

g =l 


5/ 


S' 


Substitution of [B„] and [B f ] from Eq. (A.2.7) into Eq. (A.2.8) yields 


^y Vt ) (*) = ^/^(X/)[ A 1 B w 

g = 1 


s/ 


^ ) (x) = Z^(x / ) l A “ 1 ^ 

S=1 


gr 


(A.2.6) 


(A.2.7) 


(A.2.8) 


(A.2.9) 
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A.3 C 1 Problems 

The first, second, and third derivatives of 


m 


w ( ; v) ^)-Y, p g (x j) 

g = i 


A _1 B 


W. 


gj 


m 


¥j 0) ( x ) = ^ d Pg( x j) 

g=i 


k _1 B, 


■g/ 


(A. 3.1) 


with respect to x are sought. The derivatives are found via the procedure outlined below. 
In C 1 problems, the approximations for the solution, w h (x), can be written as 

w h (x) = p T (x)a(x) , (A.3.2a) 

where p is an m >h order basis function, and a is a vector of undetermined coefficients, and 


W h (x) = ^ (Vj w) (x)Wj + ¥ { f ] (x)dj j , 

7 = 1 


(A.3. 2b) 


where y/ { p (x) and y/f Ux) are the shape functions (A.2. 9), and Wj and () j are 


fictitious nodal values. For the local weak form, the first, second, and third derivatives of 
w h (x) are needed. To evaluate these derivatives, a general procedure similar to that 
presented in section A. 1 for C° problems is used. 

A.3.1 First Derivatives 

Differentiating Eq. (A.3. 2a) with respect tox, one obtains 




dx 


dx 


dx 


(A.3 .3) 


In Eq. (A.3. 3), the second term is easy to evaluate; however, evaluation of da(x)/dx in the 
first term is not straightforward. The evaluation of this derivative requires closer 
examination and is discussed below. 
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Consider Eq. (A.2.3): 

[A] {a} = [B] {§} . (A.3.4) 

( m,m)(m ,\ ) ( m,2n)(2n,l ) 

Eq. (A.3.4) can be rewritten as 

[A] {a} =[B w ]{w}+ [B r ] {t} , (A.3.5) 

(m,m)(m, 1 ) (m,n)(n, 1 ) (m,n)(n, 1 ) 


where [B«] and [B t ] are presented in section A.2. Differentiating Eq. (A.3.5) with respect 
to x, one obtains 


(A.3.6) 

ax ax ax ax ax ax 

In Eq. (A.3.6), because the fictitious nodal values {w} and {t} are not functions ofx, the 
terms [B w ]- (d{\v}/dx) and [B,]- (d{t}/dx) vanish. Rearranging Eq. (A.3.6), 


[A] 


dx dx dx dx 


This leads to 

dx dx dx dx 

The vector {a} can be evaluated using Eq. (A.3.5), 


{a} = [A] '[BJ{w} + [A] \B 


Substituting Eq. (A.3.9) into Eq. (A.3.8), 


d{ aj 
dx 


= [A] 


l 4 Bw] {w} + [A] 


dx 


dx 


-[A] ^[A] *[B w ]{w} 

dx 


[Ar'^iAr'iBjiii 

dx 


(A.3.7) 


(A.3.8) 


(A.3.9) 


(A.3.10) 
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So, substituting Eqs. (A.3.10 and A.3.9) into Eq. (A.3.3), 


dw (x) T, , 

— 3 = P (x) 

ax 


[A] 


dx dx 


w 


-[A] 1 ^[A] 1 [B w ]{w}-[A] ^[A] \b ,]{t} 


+ 


T 

dp (x) 


dx L 


dx 


[A] ] [B w ]{w} + [A] 


dx 


(A.3.11) 


Now consider the statement of w (x) in Eq. (A.3.2b). Differentiating with respect to x, 


dw (x) 
dx 


-z 

7=1 

n 

-z 

7=1 


dyf j (x) dWj {w) - dy/j (x) dOj 

w ; — ^ +-r-¥/ (x) + 0 , — ^ + —r¥j (x) 

J J dx dx J 


7 


q?x dx 


d\f/ (x) - d y/' (x) 

w< — \-9 j 

7 dx 7 dx 

J 


(A.3.12) 


Equating the two expressions for the derivative of w(x), i.e., 
dw ,h (x) 


dx 


dw' 1 (x) 


Eq. A.3.11 


dx 


Eq. A.3.12 


(A.3.13) 


leads to 


p 1 (x){W; 


[A] l d[B w ] |-^j l </[A] r A 1 l 


dx 


dx 


[A] [B w ] 


T 

+ P (x){w 
,, ( 




dx 


dx 


+ ^(wi[ A r'[Bj 

dx 


+ Vk { i )[A] - 1 [B,] 


dx 


■z 

7=1 


dy/ ( - v \x) - di// i0) (x) ' 

Wj — + 6: — 1 

7 dx 7 dx 


V 


(A.3.14) 


Comparing the coefficients of w f and 0 j on both sides of Eq. A.3.14 gives 
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dip 


(w) m 

^=Y 

dx dx 

g = 1 


dp 


g f [A] ] [B W ] 


gj 


+ 


Pg 


[A] 


1 [A] 1 ^[A] r A 1 


dx 


dx 


[A]” [B w ] 


J gJ 


(A.3.15a) 


and 


dy/f ^ 

^=2- 


dx 


g=l 

+ p g 




gj 




dx 


dx 


J gj 


(A. 3. 15b) 


A.3.2 Second Derivatives 


Differentiating Eq. (A.3.3) with respect to x, one obtains 


0 7 O t O T 

d w (x) T d a(x) , dp (x) c/a(x) d p (x) , . 
o — = P (x) T- + 2 — + — a(x). 


dx 


dx 


dx dx 


dx 


(A.3.16) 


In Eq. (A.3.16), the last term is easy to evaluate. The term da(x)/dx was found in section 
A.3.1. The evaluation of c/ 2 a(x)/c/x 2 in the first tenn requires closer examination and is 
discussed below. 

Consider Eq. (A.3.7): 


[A] 


d{ a} _ 4B W ] 


dx 


dx 




dx 


d[ A] 

dx 


(A.3.1 7) 


Differentiating Eq. (A.3. 17) with respect to x, one obtains 
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(A.3.18) 


rA 1 /{a} , d[A] 4a}_/[B w ] ( ^ . /[ B t ] (2 
L A J — + — : 7 — 7 — ! w i+ 7 — I 1 

dx LX LX dx dx 


d 2 [ A] 

dx 2 


{a}- 


d[ A] t/{a} 
dx dx 


This leads to 


= 1 A ]' 1 ^%!(w} +1A]-' 


dx 


dx 


dx 


_2[Ar 1 4A]^}_[A]- 1 ^l { a } 


(A.3.19) 


dx dx dx 

Substitution of the expressions for {a} and d{a}/dx from Eqs. (A.3.9 and A.3.10) into Eq. 
(A.3.19) yields 

= [A ]- 1 L [ « vd { ^ } + [Af 1 dJM^} 

dx dx dx 


_2[A] 1 ^l| [A] l ^M{*} + [A] 


dx 


dx 


dx 


(A.3.20) 


-[A] l ^-[A] '[B w ]{w}-[A] 1 ^l[ A r 1 [B ? ]{t} 


dx 


dx 


[A] ’ 71 L1LL ) [A] -[BJjwJ+tA] 


n-1 d [A] I 


,-l r 


dx 


So, substituting Eqs. (A.3.20, A.3.10, and A.3.9) into Eq. (A.3.16), 
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d~w (x) T ( , 

— = p oo 

dx~ 


d~[ BJ 


[A ]' 1 + [A ]' 1 


i « [B r ] 


t/x “ 


dx~ 


_ 2 [A r'7!4lj [ A ]->^£5»l{*} +(A r>^!2iI{i} 

dx [ ax ax 

.[ A ]-i^l[Ar 1 [B»,]{w}-[A]- , ^l[Ar , [B,]{t} 


dx 


dx 


[ A]- 1 ^_M{[Ar 1 [B w ]{w} + [Ar 1 [B r ]{t}} 


dx 


+ 2 


T 

t/p (x) 


dx 


[A] 


[A] 


-1 ^[ B „ ] r . r a i-l ^[B/] 


dx 


w} + [A] 


dx 


'{t} 


- 1 ^ [A] [A ]" 1 [B w ] {w} - [A ]" 1 ^^.[A ]- 1 [B r ] {t} 


dx 


dx 


2 T 

d — l- [ [A ]' 1 [B w ] {w} + [A ]" 1 [B, ] {t} 


dx 


(A.3.21) 


Now consider the statement of dw h (x)/dx in Eq. (A.3.12). Differentiating with respect to 


d w (x) 
dx 2 


n 

-I 

7=1 


” ^ d 2 y/f ] (x) ; d 2 y/ { f\x)^ 

wj f2 — + 0j fy- 


V 


dx 


dx 


J 


(A.3.22) 


Equating the two expressions for the second derivative of w (x), i.e., 


,2 h 

d w (x) 

dx 2 


Eq. A.3.21 


,2 h 

d w (x) 

dx 2 


Eq. A.3.22 


(A.3.23) 


leads to 
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{t} [A] ‘[B,] 





[y ]p i- [ ,M a ]p i- [v]/ 7 > 



o\ 



[A] ‘[B,] 




A.3.3 Third Derivatives 


Differentiating Eq. (A.3. 16) with respect tox, one obtains 


d^w\x) T d\(x) t/p T (x) c/ 2 a(x) 
^ = p (x) 5 +3- F 


dx 


dx 


dx 


dx 


, 3 d~P ‘ (- v > da(x) | v) 


dx 


dx 


dx 


(A.3.26) 


In Eq. (A.3.26), the last term is easy to evaluate. The term da(x)/dx was found in section 
A.3.1, and the term d 2 a(x)/dx 2 was found in section A. 3.2. As in sections A.3.1 and 
A.3.2, the evaluation of d 3 a(x)/c/x 3 in the requires closer examination and is discussed 
below. 


Consider Eq. (A.3 . 1 8). Differentiating with respect to x, one obtains 


rA1 c/ 3 {a} _ r/[A]c/ 2 {a] _ c/[A]c/ 2 ja} i d 2 [A]d{a] 
[AJ - l ■ ^ I ~ H ■ 


dx 


dx dx 2 dx c / x 2 j x 2 dx 


dx dx 

d 2 [A]d{a} </ 3 [Al , d[A]d 2 {a} </ 2 [A]</{a} 

> Q ' 


dx 


2 dx 


dx 


3 18 ' dx 


dx dx 


2 dx 


(A.3.27) 


This leads to 


/{a} 

dx 


n -l d [B w ] 


= [A] ,2 -^{w} + [A] 


-1 d~ [B ? ] 


dx 


dx 


_ 3[Ar 'ii^A^_ 3[Ar >miW_T^i !a) . 

d x dx dx dx c / x 

Substitution of Eqs. (A.3. 9, A.3. 10, A.3. 20, and A.3 .28) into Eq. (A.3.26) yields 


(A.3.28) 
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[vl + WpwnTi-Cv] m 



Now consider the statement of d^w h (x)ldx' in Eq. (A.3.22). Differentiating with respect 
to x, 



and comparing the coefficients of Wj and 0 , gives 


(A.3.30) 


(A.3.31) 
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o 



(A. 3. 32a) 




([A] '[Bj) . (A.3.32b) 



Appendix B: Conditioning of Matrices 

In this appendix, the conditioning of matrices is discussed. Properties of ill- 
conditioned matrices are presented, followed by the definition of the condition number 
and an application to an example problem. 

B.l Ill-Conditioning 

An ill-conditioned matrix is one that is nearly singular, i.e., the matrix has rows 
that are almost scalar multiples of each other. Singular matrices cannot be inverted, and 
thus the inversion of ill-conditioned (nearly singular) matrices yields poor results. 

Consider a system of equations, 

[D]{R} = {P}, (B.l .1) 

for which a solution is sought. The solution is found by inverting the [D] matrix, 

{R} = [DfV}. (B.l. 2) 

In numerical applications of equation solving, an accurate computation of [D ]" 1 depends 
on the accuracy to which the components of [D] are stored, i.e., the number of significant 
digits maintained for each component of [D], The condition number of [D], cond[ D], 
provides an estimate of the number of digits lost in computing this inverse. Large 
condition numbers usually indicate ill conditioning. The question arises, “How large is 
large in terms of condition numbers?” To answer this question, the method by which 
condition numbers are calculated and an example is presented below. 
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B.2 Conditioning Numbers 

Consider the system described by (Cook et. ah, 2002) 

d\ -d\ 

— d\ d\ + ^2 

If d\ » d 2 , the second row of [D] is essentially the negative of the first row. Thus the 
matrix [D] is nearly singular and thus ill-conditioned. The conditioning number of a 
matrix [D] may be defined as (Cook et al., 2002) 

cond[V>] = ^L (B.2. 2) 

Anin 



(B.2. 1 ) 


where Anmx and 2, nm are the largest and smallest eigenvalues of [D], The eigenvalues of 
[D] are computed from 

|D-2I| = 0 (B.2. 3) 

where I is the identity matrix. In numerical computations, truncation and round-off 
errors result in the existence of errors 8 D and c>R, related to [D] and [R] by (Bathe, 
1996) 

£R = -D _1 8 D R (B.2. 4) 

Taking norms, Eq. (B.2.4) becomes 
||<?R|| \\8 Dll 

< cond[D] V/ . (B.2. 5) 

R D 


To evaluate these errors, assume t-digit precision in the computer, and .v-digit precision in 
the solution. Then, 


£D -t 

■L-J = 10 

D 


(B.2. 6a) 


and 
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(B ,2.6b) 


M 

IIrII 


= 10 


—s 


Substitution of Eqs. (B.2.6) into (B.2.5) yields an estimate of the number of accurate 
digits maintained in the solution: 

s > t-logio[co«r/[D]]. (B.2.7) 

In the system of Eq. (B.2.1), consider d\ = 2 and d 2 = 1 (cl\ > d 2 ): 


[D] = 




(B.2.8) 


To compute the eigenvalues. 


" 2 -2 


~ A 

0 “ 


1 

1 

1 

K) 
1 

-2 3 _ 


0 

A_ 


l 

1 

bo 

1 


and, according to Eq. (B.2.3), 

2-/1 -2 

, =0 

-2 3-/1 

(2 - A)(3 - A) - 4 = 0 . 

The eigenvalues are therefore 
ii = 4.56155 
A 2 = 0.43845 

and the condition number is calculated as 


(B.2.9) 


(B.2.10) 


(B.2.1 1) 


cond[D]= 4 56155 = 10.4038. (B.2.12) 

0.43845 

Assuming a double precision computer is used, t = 14, and the number of accurate digits 
maintained in the solution can then be computed as 

5 > 14 -log! 0 [l 0.4038] = 13 . (B.2.13) 
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Now consider the case d\ » eh, for example, d\ = 2 and cfe = lxl 0 °: 


[D] = 


2 

-2 


-2 

2.000001 


From Eq. (B.2.3), 


2-2 -2 
-2 2 . 000001-1 


(B.2.14) 


(B .2.15) 


and the eigenvalues are 

21 = 4.0000005 (B .2.16) 

2 2 = 0.0000005 . 

The condition number is therefore calculated as 
Jrrv , 4.0000005 0 AAA AA t 

cond[D] = = 8000001. B.2.17 

0.0000005 

Using the same double precision computer, t = 14, the number of accurate digits 
maintained in this solution is computed as 

s>14-log 10 [8000001]= 7. ( B .2.18) 

In the entry 2.000001 of Eq. (B.2.13), the “1” that keeps the matrix from becoming 
singular is in the seventh significant digit location. Because only seven digits are 
maintained during subsequent computations (see Eq. B.2.17), the inverse of [D] in Eq. 
(B.2.14) will be very inaccurate. Thus, the conditioning number of a matrix is a good 
indicator of how well the matrix is conditioned. 
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